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NOMENCLATURE 


heat capacity 

constants in heater output 
function 

enthalpy 

heat transfer coefficient 

thermal conductivity 

latent heat of fusion for ice 

layer thickness 

rate of heat production per 
unit volume 

rate of heat production per 
unit area 

temperature 

time variable 

heater time on and time off 


(Btu/lb-°F) 


(Btu/ft 3 ) 
(Btu/ft 2 ~hr-°F) 
(Btu/ft-hr-°F) 
(Btu/lb) 

(ft) 

(Btu/hr-ft 3 ) 

(Btu/hr-ft 2 or 
Watts/in 2 ) 

(°P) 

(hr) 

(hr) 


dependent variable, H or T 

fraction of Diodal volume 
which is ice 

space coordinate in one- (ft) 

dimension 

fraction of nodal volume 
which is liquid water 

position of solid-liquid (ft) 

interface 


Greek letters: 


Ot thermal diffusivity 

At time step 

ax grid spacing 

P density 

U) over-relaxation parameter 


(ft 2 /hr) 

(hr) 

(ft) 

(lb/ft 3 ) 


Subscripts : 

al inner ambient boundary 

a2 outer ambient boundary 

i layer in composite blade 

II outer layer of composite blade 

(abrasion shield) 

j grid point 

1 liquid (water) 

Imp liquid at the melting point 

mp melting point 

s solid (ice) 

smp solid at the melting point 

w ice-water layer 

Superscripts : 

• point heat source 

o evaluated at the previous time 

level 

a evaluated halfway between the 

previous and present time level 

(old) value from previous iteration 

(new) value from current iteration 


v» 


IV 
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I . INTRODUCTION 


The formation of ice on aircraft components poses a 
problem of considerable significance. For the aircraft to 
perform safely and efficiently at near or below freezing 
temperatures, this ice must be removed. Both anti-icing and 
de-icing systems are used for this purpose. An anti-icing 
system prevents the formation of ice, whereas a de-icing 
system periodically removes the ice that has formed. This 
investigation deals with electrothermal de-icing as applied 
to ice removal from propeller and helicopter rotor blades. 

As such, this is a continuation and extension of the work 
done by G. Baliga [1] at the University of Toledo. 

A de-icer works by destroying the adhesion between the 
ice and the composite blade surface, thus allowing aero- 
dynamic or centrifugal forces to sweep away the ice. This 
is accomplished in an electrothermal de-icer pad by means of 
a resistance heater which raises the temperature of the com- 
posite blade surface above the melting point of ice. Since 
only a thin layer of ice need be melted to destroy adhesion, 
the energy requirements are significantly less than those of 
other systems. Baliga [1] has reviewed the advantages and 
pitfalls of other anti-icing and de-icing systems. 

A section of an electrothermal de-icer pad embedded in 
an aircraft blade is shown in Figure la. It is a composite 
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body consisting of five layers. The center layer is the 
heater, which is separated from the substrate and the 
abrasion shield by insulating layers. The bonding between 
these layers is suspected bo be less than perfect, and small 
air pockets between layers may exist. In operation, the 
heater is turned on periodically to remove ice that has 
formed on the abrasion shield surface. 

Typically, the heater is a woven mat of wires and glass 
fibers or multiple strips of resistance ribbon. Woven mats 
may have thicknesses as great as 0.020", whereas ribbons 
have thicknesses between 0.001" and 0.005". Individual 
heating elements are between 0.5" and 1.0" wide. Stalla- 
brass [ 2 ] has pointed out that gaps which exist between 
these heating elements can reduce the effectiveness of the 
de-icer pad, causing non-uniform of the ice. The 

gap width is roughly 0.080" for woven mats and 0.040" for 
metal ribbons. 

The two layers adjacent to the heater provide electrical 
insulation. In order to direct most of the heat outward, the 
outer layer should have a much higher thermal conductivity 
than the inner layer. This is generally not possible since 
good electrical insulators are also poor conductors of heat. 

To compensate for this effect, it is necessary to use a much 
greater thickness for the inner layer. A ratio of thicknesses 
of at least 2:1 has been recommended. Resin- impregnated 
cloth is commonly used for both layers. Electrical insulation 
requirements necessitate that the outer layer of cloth have a 
thickness between 0.010" and 0.020". 
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The purpose of the abrasion shield is to protect the 
de-icer from the environment and also to cut down drag on 
the composite blade surface. For these reasons , stainless 
steel is normally used for the abrasion shield. The rela- 
tively high thermal conductivity of stainless steel enables 
heat to be conducted laterally across the blade. This can 
be beneficial to melting the ice above the heater gaps. 
Thicknesses for the abrasion shield range from 0.010*' to 
0 . 020 ". 

A wide range of materials is used as substrates 
depending on the particular application. An aluminum alloy 
is considered in this study. Baliga [1] and Stallabrass [2] 
have examined the effects of different materials and thick- 
nesses on de-icer performance. 

Due to the large number of parameters that affect the 
rate of heat transfer in the composite blade, it is not 
surprising that many proposed de-iaer designs fail to 
achieve the level of performance expected. This complexity 
leads naturally to the use of numerical methods along with 
the digital computer to evaluate de-icer performance. In 
this study, de-icer performance is measured by the time 
required to melt the ice at the ice-abrasion shield inter- 
face (or a finite thickness of ice) starting from various 
initial temperatures . The model constructed considers one- 
dimensional, unsteady-state heat transfer in a composite 
body. A wide range of parameters are available to completely 
specify the de-icer design. The phase change in the ice 
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layer is accounted for by the Enthalpy method. This method 
and the numerical methods employed in the model are reviewed 
in the next section. The complete numerical formulation of 
the problem appears in Section III. 
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II. LITERATURE REVIEW 

There have been several recent studies concerned with 
the performance of electrothermal de-icer pads. Of these, 
only the investigations of Baliga [1] and Stallabrass [2] 
have considered the effects of the phase change in the ice 
layer. Gent and Cansdale [3], while not considering the 
phase change in their simulation, do present temperature 
profiles from experimental de-icer pads. The de-icer pad 
model used in the present study takes into account tire phase 
change, and also contains significant improvements over the 
models used in the studies mentioned above. All of these 
models have been one-dimensional except for that of Stalla- 
brass, who also developed a two-dimensional model. The 
analytical and numerical methods used in the present and 
previous studies are outlined below. 

A. ANALYTICAL METHODS 

A variety of analytical techniques is available to 
solve transient heat conduction problems in composite 
bodies, the most common being the Laplace transformation. 

However, most of these techniques are too complicated to 
apply when the body contains more than two layers. An 
exception is a method proposed by Campbell [4], where the 
analogy between one-dimensional heat conduction and the flow 
of electricity along a transmission line is used to calculate 

' ft 

... .... . 2 . 
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the temperature at any point within the composite body. 
Stallabrass [2] used Campbell's method to check the accuracy 
of his numerical technique. All analytical methods, however, 
have the disadvantage that an excessive amount of calcula- 
tions must be done for each temperature desired. In 
addition, they cannot be used when the pi, vise change in the 
ice layer is considered. 

B . NUMERICAL METHODS 

All of the recent models proposed for an electrothermal 
de-icer pad have used finite difference methods. In these 
methods, the differential equation governing the heat trans- 
fer in the composite body is replaced by a system of 
difference equations. This transforms the continuous time 
and space domain of the problem into a discrete grid of nodal 
points. The difference equations can then be solved alge- 
braically to determine the temperature at all nodal points 
at any time step. This is a definite advantage over analyti- 
cal methods. Finite differencing is an approximate technique, 
and its accuracy depends upon which of the several finite 
difference schemes is used along with the grid spacing (ax 
and At) chosen. The accuracy is measured by the order of the 
truncation error for both the time and space derivatives. 

For some of these schemes, a restriction also exists on the 
size of ax or At that will ensure convergence and stability 
of the solution. A finite difference representation of a 
de-icer pad appears in Figure lb. 
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Both Stallabrass (2 ] and Gent and CansdaXe (3) used the 
explicit forward finite difference scheme. In this scheme, 
the temperature at a node can be calculated directly from 
the nodal temperatures at the previous time step. The 
truncation error is first order in time and second order in 
space. The convergence and stability criteria for forward 
differencing is: 

(X At/(Ax) 2 $l/2 where CX is the thermal diffusivity of the 
layer in the composite body. For the de-icer problem, this 
requires a time step of 0.001 sec. or smaller to be used. 

The excessive number of calculations needed because of this 
small time step can cause an accumulation of truncation and 
round-off error. 

In Baliga's work [1] and in this study, the Crank- 
Nicolson implicit finite difference scheme is used. This 
method is unconditionally stable and no restrictions are 
placed on the size of At and ajc. In addition, the truncation 
error is second order for both time and space. This allows a 
time step of 0.1 sec. to be used, thus reducing the total 
number of calculations . The only drawback of this method is 
that the temperature at any grid point can no longer be 
explicitly calculated. The system of equations which results 
must be inverted or else solved iteratively in order to 
obtain the temperature distribution at any time step. Baliga 
used the method of Thomas to invert the tridiagonal system of 
equations. The method employed for the phase change in the 
present study dictates that Gauss-Seidel iteration be used. 
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This requires more calculations to be done but reduces 
round-off error, 

C. METHODS FOR HANDLING THE PHASE CHANGE 

In the past few years, there has been a significant 
increase in the number of articles appearing in the litera- 
ture that deal with phase change and related moving boundary 
problems. These types of problems are sometimes referred to 
as Stefan problems. Due to the nonlinear boundary condition 
caused by the movement of the solid-liquid interface, these 
problems are relatively difficult to solve. Analytical solu- 
tions are only available for simple problems and many numerical 
techniques have been proposed. An extensive review of most of 
the analytical and numerical techniques that have been used 

in Reference 5, Many of these methods use predictor- 
corrector techniques, where the phase change interface loca- 
tion is assumed, and subsequent iterative calculations correct 
this position. This requires an excessive amount of calcula- 
tion. The added complexity of the heat transfer occurring in 
the rest of the composite body makes these methods impracti- 
cal for the de-icer problem. For this reason, methods which 
do not require trial and error calculations to determine the 
interface location have been used. 

Stallabrass [2] accounted for the phase change by holding 
a node at the melting point until enough energy had been 
transferred to completely melt the nodal volume. Baliga [ll 
approximated the latent heat effect with a large change in 
heat capacity over a small temperature interval around the 
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melting point. The thermal conductivity was also allowed 
to vary linearly over the interval. This technique was 
proposed by Bonacina, et. al. [6] . Both methods are very 
similar to the Enthalpy method, but lack the formalism 
which makes this method easy to apply numerically. The 
Enthalpy method, which is also called the method of weak 
solution, is used in this investigation. 

In the Enthalpy method, the governing equation for 
conservation of energy is formulated in terms of two de- 
pendent variables, enthalpy and temperature. The moving 
boundary condition and predictions of the phase change inter- 
face location are not needed. After the enthalpy at a node 
is calculated, the known enthalpy-temperature relationship 
for water can be used to determine the nodal temperature. 

The equivalence of this method to the moving boundary formu- 
lation was proven by Atthey [7] . 

Most of the applications of the Enthalpy method have been 
formulated using the forward finite difference scheme. In 
this study, the Crank -Nicolson scheme is used, and the system 
of equations which results is solved by Gauss-Seidel itera- 
tion. Voller and Cross [8,9] in two recent articles have 
pointed out that the Enthalpy method yields unrealistic 
results since a node remains at the melting point for a 
finite period of time. This leads to the prediction of 
temperatures which oscillate around their true values. The 
same phenomenon also occurs with the methods of Stallabrass 
and Baliga. By reinterpreting the Enthalpy method, Voller 
and Cross have derived a criterion for determining the points 
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of correspondence between the true and oscillating curves. 
This enables accurate temperature profiles to be obtained. 
The criterion is given in the "Discussion of Results" 
section. 
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III, NUMERICAL FORMULATION 

A. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 

Tlie following assumptions were made in the formulation 
of a one-dimensional, unsteady-state, mathematical model for 
heat transfer in a composite aircraft blade on which an ice 
layer has formed : 

1) The physical properties of the materials composing 
each layer of the composite blade are independent 
of temperature; 

2) Lateral heat transfer in the layers can be 
neglected, so that only a one-dimensional model 
need be constructed; 

3) The ambient temperature and all heat transfer 
coefficients are constant; 

4) The ice layer thickness is constant; 

5) The effect of the volume contraction of the ice 
as it melts can be neglected; and 

6) The ice is "pure", so that the latent heat is 
released isothermally at the melting point. 

1 . Composite Aircraft Blade 

With the above assumptions, the governing differential 
equation for each layer of the composite aircraft blade is : 

Pi ^pi = k i +Q i = l,...,n (1) 
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where 


T i 


Q. 

X 



temperature in layer i 

rate of heat production per unit volume in 
layer i 

density of the ith layer 

heat capacity per unit mass of the ith layer 
thermal conductivity of the ith layer 


x = space coordinate 
t = time variable 

II = number of layers in the blade 


A composite blade containing a finite thickness heater is 
characterized by: 


i = 1 

t 

substrate 


© 

It 

O 

1 = 2 

» 

lower or inner 

insulation 

Q = 0 

2 

i = 3 

i 

heater 


Q = Q (t) 

3 3 

tl 

t 

upper or outer 

insulation 

© 

II 

O 

i = 

II 

=5, abrasion 

shield 

© 

ji 

II 

o 


A variety of different boundary conditions is considered 
with equation (1). These are: 

(i) For perfect contact between layers, the temperature and 
heat flux are continuous at the layer interfaces. This 


leads to the boundary conditions: 


T i I = T i+1 


(3a) 


i=l, . . • ,11-1 


-k 


= -k. , . £Si±l 
I 1+1 d)x 


(3b) 


where "I" denotes an inter facts 


(ii) In reality, there mty exist a resistance to heat trans- 
fer across the layer interfaces due to the small layer 
of adhesive used to hold adjacent layers together and 
also to small air gaps caused by poor contact. The 
boundary conditions for this case are: 




= h. (T. 
I i i 


- T 

I i+1 


li> = 


(4a, b) 


i—l, . « « , II— 1 


where h^ is the heat transfer coefficient across an 
interface . 


(iii)If the heater can be treated as a point heat source 

(zero thickness), an alternate equation is used for the 
interfaces between layers, which is: 



where is the rate of heat production per unit area. 
Equation (3a) still applies at an interface. 

A blade with a point heat source is characterized by: 


i = 

1, 

substrate 

Q x = o. 

= o 

i = 

2, 

inner insulation 

Q a = o. 

q ' =q'(t) 
* 2 (6) 

i = 

3, 

outer insulation 

Q = 0, 
3 

q ' = 0 

3 

i = 

II 

= 4, abrasion shield 

© 

il 

O 

V 

o 

II 

<1 


4 4 


(iv) Convective heat transfer occurs at the inner boundary 
of the composite blade and also at the outer boundary 
if the ice layer is not present. For the inner 
boundary : 
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where "1" denotes the inner ambient boundary, h a ^ is the 
convective heat transfer coefficient at the boundary and 
T al is the ambient temperature. Since the air within 
the blade is stagnant, h_^ is small. 

(v) For the outer boundary: 


J-55T 


= h 


a2 


(T i 2 “ 


a2 


i : =II 


( 8 ) 


where "2” denotes the outer ambient boundary, h a2 is the 
convective heat transfer coefficient at the boundary and 
T a2 is the ambient temperature. The quantity h a2 is 
very large due to the dynamic forces acting on the out- 
side of the blade. 

Besides the above, constant temperature boundary condi- 
tions can be specified for the inner and outer surfaces of the 
composite blade. The initial temperature distribution in the 
composite blade can be constant or a function of position. 

2 . Heat Source 

The total output of the heater is the same regardless of 

whether it is treated as being of finite or zero thickness. 

' 

Thus, the total rate of heat production per unit area is: 

q i (t) = ljO^t) = q-jL^t) (9) 

where 1^ is the thickness of the heater. A wide range of 
different heater outputs can be specified. These include: 
outputs that are constant, linear or sinusoidal with time, 
and also outputs that can be periodically turned on and off: 
ramps, square waves, etc. The general expression for these 
functions is: 
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q i (t * 


’ 

, c,t + c + c cos (c t 

1 2 4 

0 


q ± (t+P) = q i (t) 


+ c 5 ), 0<tst on 

' t on <t <P 

Pst + t, 

on off 

( 10 ) 


t > P 


where C., C , C , C and C are constants, t and t are 
JL 2 3 4 5 on Off 

the times the heater is on and off, respectively, and P is 
the period of the output. 

3 . Ice Layer 

The classical formulation for the ice layer subject to 
assumptions (4), (5) and (6) is: 

P/ps ^ = * s I^Sa. x>y Ula) 

s * <9t s 




x < y 


(lib) 


along with the moving boundary condition: 


To = T, = T 
s 1 mp 


x =y 


(12a) 


c^Ts 

53T 


-h. 


d)Tx 


= PL Jjif 
y dt 


(12b) 


where 


T s = temperature within the solid 

T^ = temperature within the liquid 
T = melting point 

P g , ^) S ,k s = physical properties of the solid 
Pi* = physical properties of the liquid 

pL = latent heat of fusion per unit volume 
y* = position of the solid-liquid interface 
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Ab discussed in the "Literature Review", the solution of 
equations (11) and (12) requires that the interface location 
be solved for explicitly. To avoid this difficult procedure, 
the Enthalpy method is applied. The governing differential 
equation for the Enthalpy method is: 






w _ <3 


35 


( k 


(13) 


where 


«w = 


enthalpy per unit volume within the 
ice-water layer 


T = temperature within the ice-water layer 


w 
k w = 


thermal conductivity within the 
ice -water layer 


Thus, the enthalpy within both phases is found using only 

one equation. The known H vs. T relationship is used 

w w 

to determine T ? this relationship is: 

W 


/)$ T 
r'a ps w 


, T < T 
9 w mp 


H w = 


CL'S , ( T - T )+£.((? T + L) , T > T 
n. pi w mp ; “l' ps mp ' * w i 


(14) 


ps mp 


mp 


where L is the latent heat of fusion per unit mass. It has 

been shown elsewhere [7] that the formulation above is 

equivalent to the moving boundary formulation, equations 

(11) and (12). For numerical solutions, it is easier to 

work with T as a function of H . Inversion of (14) gives: 
w w 
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( VPs e , 


ps 


' H W < H smp 


T = ( 
w 


mp 


H 


smp 


< H < H 


w 


Imp 


(15) 


1 < H w- H lmp>/Pi e pl + T mp * H > H 


Imp 


with 


H 


smp 


= 0$ T 
"s ps mp 


H 


where H 


Imp 
and H 


= P (S' T ■* L) 
“1 ps mp ' 


ps mp 

are the enthalpy of the solid and the 


smp Imp 

liquid at the melting point, respectively. Also note that 
in equation (13), the thermal conductivity is now a function 
of position. 

Boundary conditions must also be specified for the ice- 
water layer at the interfaces with the abrasion shield and 
the atmosphere. Perfect contact between the layer and the 
abrasion shield is assumed, so that equations (3a,b) apply 
with i+l=w, which ares 


T . = T 

i I w |l 


1 — XX 


(16a) 


-k. <$Ti 


_ -*w c9T w 

^r 


(16b) 


Equation (8) holds for the outer boundary of the ice-water 
layer with i=w, which is : 


2 = h a2 ( t „I 2 - T a2 > < 17 » 

After a thin layer of ice has melted, the layer can be 
shed by the dynamic forces acting on the composite blade, 
and equation (8) applies at the outer boundary. 
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B. CRANK-NICOLSON FINITE DIFFERENCE FORMULATION 


For numerical solutions, the above differential 
equations are replaced by their finite difference analogs. 
In this study, the Crank -Nicolson finite difference scheme 
is used. 

Truncated series expansions are used in the Crank- 
Nicolson scheme to approximate the partial derivatives 
appearing in the governing equations. Letting U stand for 
the dependent variable, either T or H, truncated Taylor's 
series expansions for the partiai derivatives and 0 
are: 


and 


do 

3S 


d 2 U 

dx2 


= u i+l - pi-1 + oUx) 2 

j 2 ( a x) 


_ U i+ i - ?U. 

g- (Txl* 


+ Uj-1 


+ OUx) 


(18) 

(19) 


where the subscripts j-1, j and j+1 denote adjacent nodal 
values. The grid spacing, ax, is constant within a layer, 
but may vary between different layers. Equation (13) requires 
the ejqpansion for (k ) f which is : 



fcj+y^Uj+l ~ u j) - ^j-i/ x ( u j~ u j-l) 
(ax) 2 


+ 0(ax) 2 

( 20 ) 


where k . , and k , are average values of k between nodes 
3+W 3-i/x 

j+1 and j, and nodes j and j-1, respectively. The truncation 

error for these approximations of the partial derivatives is 

second order. The second order finite difference analog for 

the time derivative, s£H, is: 

3t 
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A u *i “ u j° 

= - 2 __ i + 0( At) 2 (21) 

j At 

where the superscript 0 denotes the value at the previous 
time level and the superscript t denotes the value halfway 
between the previous and present time levels. The time step 
At can be changed as the calculations progress in time. 

In the Crarik-Nicolson scheme, the governing differential 
equations and boundary conditions are approximated at a point 
halfway between the known and unknown time levels. 'The 
approximation for the time derivative is given by equation 
(21) . The analogs for the space derivatives given in equa- 
tions (18, 19 and 20), however, cannot be used since they 
would require the. evaluation of the dependent variable at 
the half time level. To overcome this difficultly, the 
Crank-Nicolson scheme employs the following approximations 
for these derivatives : 





Vi ‘ 2u j + Vi + U 1°+1 


2 (ax) 2 


- 2U.° + 




u ° 

_Jhi 


+ O(Ax) 2 (23) 
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(24) 

”j°-l> 


2 (Ax) 


+ O(ax) 


In addition, the approximation for U. evaluated at the half 

J 

time level is : 


Uj = ~ (U,. + U°) + 0 ( At ) 2 


(25) 


The Crank -Nicolson finite difference equations are obtained 
by substituting the above into the governing differential 
equations and boundary conditions of Part A. 


C. FINITE DIFFERENCE EQUATIONS FOR THE COMPOSITE AIRCRAFT 
BLADE 

Substitution of the analogs (21) and (23) into equation 
(1) yields: 

£$• Tj a — T j . -i 


p.e . j i = 

' i pi 

At 

k. ' 2 T i.i + + Ti ° 3 +1 ' 2 Tj °>3 * Ti ° 3- 

2(ix.) J 


(26a) 


+ Q i 


• • • f XX 


Solving for the unknown temperature at node j gives s 


T. . = (T. + T. . , + T.° ... + 2(M.-1) T.° . + 

1,3 V 1,3+1 i, D“*l i, 3+1 i ’ ~ 


i»3 


T i,j-i + 2 s i)/ 2 (M i +1) 


( 26 b) 
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where M i = U^) 2 /& i At 

s i = Q i ( 4 x i )2/k i 

«L = V^pi 


The quantity is the thermal diffusivity of the ith layer. 
The source term, S . , is a function of time and so is evaluated 

i»» 

at the half time step. Equation (26) is valid for all grid 
points within each of the layers, 1 through II. An alter- 
nate expression is developed later for the ice layer. At the 
interfaces between layers and at the inner and outer surfaces 
of the composite blade, the boundary conditions must be 
finite differenced. The finite difference analogs of condi- 
tions (i) through (v) of Part A are given below: 


1. Perfect Contact Interface - B.C, (i) 

For this case, let j be the inter facial node between the 
layers i and i+1 as shown in Figure 2a. Finite differencing 
equation (3) with the aid of analog (22) gives: 


T. . = T . , . . , T.° . = T° 

1,3 1+1,3 ' 1,3 l+l. 


(27a) 


1*— 1 , c . . , II 1 ** 


and 


T - T + T ° — T ° 

i,3+l i,j-l i, j+i i,j-l _ 

4axj 


<3 

m —fp + mO ^ m O 

-k. +1 i+1,3+1 i+l, j~l A i+l,J+l i i+l,j-l 

4Ax i+l 


(2 7b) 
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The nodal temperatures T ljj+l , TjP j+1 , T. +1 ^ and TP +1#j-1 

are fictitious values and must be eliminated from the finite 
difference expression. This is done by the application of 
equation (26) to node j for both layers i and i+1. This 
yields : 


<r» 4. m o 

i, j+1 T i,j+1 


2( M i+ l) T 1#j - T.^^ - 2(M.-1) TO, - T f - 2S. (28a) 


T.,. . . + T.° 

i+1, 3-1 x+1,.3-1 


2 (M i+l +1)T i+l, j~ T i+l, j+1' 2 (M i+l“ L) T f+i, j“ T i°+l, j+l“ 2S i+l 


(28b) 


Equations (27a,b) and (28a,b) can be combined to eliminate 
the fictious temperatures, yielding: 


T. 4 


+ H i T i + i, j+ i + + + 

T.o . + N. T 1?+1) j+1 + S. ♦ N iSl+1 )/ 


where 


(U+M^ + ETU+M^)] 


N i = k i+1 ix i A i A>c i+1 


(29) 

i=l, . . .,11-1 


2 . Resistive Interface - B.C. (ii) 

Let j be the inter facial node for layer i and j+1 be 
the interfacial node for layer i+1. This is shown in 
Figure 2b. Substituting the analogs (22) and (25) into 
equation (4a) gives: 
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. kj T i,j+i ■ + T i,j+i " 

4Ax. 

JU 

h i(j * T i, j + T i i j ) “§ (T i+l,j+l * T i°+l,j+l ) ) 


(30) 


** 1 , # 0 • , I 1 1 

Using equation (28a), the fictitious temperatures T. ... and 

i, 3+1 

T.° + can be eliminated. Thus, the equation for T. .at this 
*i j tj> i,3 

interface is, after rearrangement: 

T. . = ( T. . - + R. . T . . ... + T.° . . - (1-M. + R. . ) T.° . 

x ,3 V li x+ 1 , 3+1 x, 3-1 x lx' A i, j 

\ / (31) 

+ R ii T i° + i, j+ i + s x )/ (1 + “i + R u> 


where 


R. . = Ax i h i 

■LX k 

x 


Similarly, for boundary condition (4b) 


T — T +T° -TO 

-k. +1 i+ltj+2 A i+1, 3 x+l,j+2 i+1, j 

4Ax i+i 


(32) 


h i& (T i,3 + T i°.3 > -K + l. j+ l * T i° + l, j+ l>) 

i=l, . . . , II-l 

Equation (28b) written for layer i+1 is used to eliminate 
T i+l,j and T i+l,j- This yields: 

T i+l,j+l = ( T i+l,j+2 + R 2i T i,j + T i+1, j+2 ' (1 '' M i+l +R 2i ) ' 


T i+l,j + l + R 2i T i,j + S i + l)/‘ 1+M i + l +R 2i> 


( 33 ) 
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where R 2i = AX i+l h i 

k i+l 


3 . Point Heat Source - B.C. (.iii ) 

The same grid as used for (i) (Figure 2a) is used for 
this interface. Applying the Crank -Nicolson derivative 
approximations, the finite difference analog of equation 
(5) is: 

-k . i.j+i i;j-l j+1 i^j-1 + a i = 

1 - ■ ■ “ ‘ 


4ax. 


T - T + TO 

-k i+i i+l,j+l i+l,j-l x i+l, j+1 


(34) 


T P 


i+1, j-1 


4Ax i+i 


i—l , • • * , X I “1 


Equations (28a, b) with S.=S. ,,=0 and T. . = T. - . and T.° .= 

1 l-rl IjJ ^*>3 

^i+1 j are Use< ^ fco e liminate the fictitious temperatures 

T. . . , , TP ... , T. , . _ and T.° . The finite difference 

i,j+l i,j+l i+ 1 , 3-1 x+ 1 , j -1 

egression for T^ ^ at the interface is : 

T i,j=( T i,j-i + N i T i + i, j+ i + + «v 1> + 

Hi (M i+ i-l)] TO + N. T.o +ljj+1 + 2 S<)/ (35) 


[(1+M i ) + N ± U+M i+1 )] 


s i = q l ix i /k i 


where 
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Like S., S'/ is evaluated at the half time step. The 

«L*- X 

similarity of this equation and equation (29) allows them 
to be combined for computer implementation. 


4 . Inner Ambient Interface -B.C. (iv) 

The grid used for the boundary between the substrate 
and the interior of the composite blade is shown in Figure 2c . 
At this surface i=l, j=l. The finite difference form for 
equation (7) is; 


T -T +T° _ T ° 

1 hi X >° l* 2 

4 Ax, 


= h al(f CT l,l + T f,l> - T al)‘ 36 > 


Temperatures T^ Q and T° q are fictitious and are eliminated 
by the same procedure as used for the other boundary condi- 
tions, this gives for T^ 

T l,l = ( T l,2 - + “al> *1.1 + T l°,2 + 2 “al T al) < 37 > 

/ (1 + M 1 + "al> 

where N , = Ax.h ./Tc. 

al 1 al 1 


5. Outer Ambient Interface - B.C. (v) 

When the ice layer is not present, the grid at the outer 
surface of the abrasion shield is that shown in Figure 2d. 

At this surface i=H, and let j be the inter facial node. 

Finite differencing equation (8) yields; 
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- T„ 


+ T, 


- T, 


-kjj IX,j^l H,j-1 ' *II,j+l _ 


4Ax 


II 


h »4 <, II.J + *n.i> - T a 2 ) 

The fictitious temperatures, T ri ^ +1 and T^. .. +1 , are 
eliminated as previously described to give: 

*11. J =( T II,j-l - (1 - *XI + N a 2> T II,j + *i # x,j-l 

+ 2 N a2 T a2 )/< 1 + «II + N a 2> 


( 38 ) 


(39) 


where N &2 = AXjjh^Aj.j 


D. FINITE DIFFERENCE EQUATIONS FOR THE ICE LAYER 


Unlike the composite blade, two equations are needed at 
each node of the ice layer; one to calculate enthalpy and 
one to calculate temperature. Substitution of analogs (21) 
and (24) into equation -(13) yields the difference equation 
to be used in the Enthalpy method,: 


H . -H ° . 
w,3 w , 3 

At 


= ( k w, j+i/,( Tw *3 +1 “ T w,j) ~ ^ j-w,( T w,3 - 


T w,j- 1 ) + j+t/ t ( T w, j +1 “ Tw <j 

k w, j -i/j T W, j ~ T w°,j -10 (AV 2 


) - 


(40) 


The equation above must be solved explicitly for the nodal 
enthalpy, H . . This requires equation (15) to be used to 

W 9 J 

relate T . to H . . Note that this leads to three sets of 
' w > J w *3 

equations; one each for the node below, at, and above the 
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melting point. Substitution of (15) and rearrangement 
yields, for the solid range: 


! w,j ( H w,j * V 2 tk w, j+1 + k w, j-./ a T w, j-1 + k w,j+i/ x * 

(T w°,j+l “ T w°, j ) - k w,j-, /x CT w,j “ T w, j-1? ^ 

[1 + £ pg ) ( k Wj j +l/i + k W,j-l /x )J 


H . < H 
w, j ** smp 


with 


T . = H ,/P 6 
w,3 vr,yr s ] 


ps 


(41b) 


for the melting range: 


H w, j = H w°,j + V 2 [k w, j+^ T w, j+1 “ ^p* " k w, j-^' 
^ T mp " T w,j-1^ + k w, j+lJ T w, j+1 “ T w,j* 


(42a) 


- T w°, 3 -i> > 


H smp < H w, j < H lmp 


with T . = T 

w,3 mp 


(42b) 


and for the liquid range: 

H , = (hO . + M/2[k ... (T ... - T + H. / 

w, j ' w, 3 w' w, 3+l/v w, 3+1 mp Imp 

P.S .) - k . , (T - H /p.e . - T . ■ ) 

1 pi w, j-i/T, mp Imp 'l pi w, 3-I 

w, j+|/^ w, j+1 w, j w, w, j 

v.j-i’M 1 + i<«w/Pi® P i )(k w.jV Vh 11 


(43a) 


H w, j ^ H lmp 
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with T w,j =(H w,j - + V 

where = At/( 4 x„) 2 


(43b) 


The algorithm used to implement these equations is presented 
later. When node j is solid at both the present and pre- 
vious time levels, equations (41a,b) can be combined. The 

resulting equation for T . is equivalent to equation (26b) . 

w, J 

This equivalence also holds when node j is liquid at both 
time levels and (43a, b) are combined. 

Equation (42) has a different form from (41) and (43). 
The node is held at the melting point, and heat entering the 
nodal volume is used only for melting. Thus, the fraction 
of the nodal volume melted can be related to the enthalphy 
of the node calculated using (42a) . Letting X^ be the 
fraction of the node which is solid and Yj be the fraction 
of the node which is liquid, an energy balance yields : 


H . = H X . + H, Y, 
w, 3 smp 3 Imp 3 


(44) 


where X. and Y. are related by: X,. + Y. =1. The movement 

3 3 * 3 3 

of the ice-water interface through the layer can be fol- 
lowed by using equation (44). 

Equation (42) can also be derived by finite differencing 
the moving boundary condition (12) and applying (44). 
Equation (12) is applied to a single node, and then finite 
differenced to give: 


T . = T 
w, 3 mp 


(45a) 
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k 


<9 T w 

w <57" 


A 

j + 


-k ^ T w 


= P L 

D- 



(45b) 


The approximations for the first derivatives with respect 
to x used in (45b) are slightly different than equation 
(22). They are: 


k d T w 
w 57“ 



k ( T 

W,j 


T w, j- 1 * + k w, j-l/i T W, j 


T 0 . ) 


(46a) 




v 

37“ 


j+ 


k ... (T 
w, 3+1/* w, 3+1 


m \ j. v O /m O 

W, j ' W, j+l/a W, j +1 


T° .) 

w,3 


(46b) 


2A \, 


The latent heat of melting per unit volume, pL, is (neglecting 
volume contraction) equal to (H^^ - H g ) ,. Using analog (21) 
for the time derivative and substituting the above finite 
difference analogs into equation (45b) yields: 


- H ) y j - y j° 


Imp smp 


At 


= kw ' 3 + W (T w, j+1 ~ T mp ) + 
2A *w 


j + Ux T w, j+1 “ T w, " k w,j-i^ T mp ” T w, j-l* 


( 47 ) 


30 


Dividing equation (47) by ix , the left hand side of the 

W 

equation becomes : 


, . y./ax - y?/Ax, 

(H_ - H ) D w 3 w 

' Imp smp rr 

At 


, Y. - Y ? 
< H lmp - H smp> -L^-1 


Wj - H s m p < 1 - V - H lmp Y j + » SmP (1 - X P (4B) 

At 


<Wj + H lmp Y -i> - H smp + H smp ~ (H smp X j + H lm P V 

At 


H . - H . 
w, 3 w , 3 


At 


In the above, equation (44) was used along with the fact that 

Yj, 

= — -J-. Substitution into equation (47) yields: 


H w,j ~ H w,j 

At 


k w, j+i/^w, j-H ~ T mp^ " k w, j-i/| T mp ~ T w, j-1* 

2 (AX w )2 


1 vo ( T o - T 0 . ) - k ° . (T°. 

K W, j+l/x w, j +1 W,j/ w,K 


m O \ 

T W, j-l' 


or 


H 


1 W, j H w, j + V 2 ^ k w, j+1/i T w, j+1 _ T mp* - k w, 


%p - T w,j-1 ! + k w, j+1/t (T w°,j + l - ^ 




The final result is equation <42a) . From the above analysis, 
it is apparent that the Enthalpy method is a much easier and 
more direct method to derive the melting point equation. 

This analysis, however, does show the equivalence of the 
classical and Enthalpy method formulations. 

As noted earlier, in the Enthalpy method formulation, 
k w is a function of position. More specifically, it is a 
function of the liquid-solid interface position. Equations 
(41) through (43) require average values for k between 
adjacent nodes. A volume average is used to ensure that the 
correct values are obtained when X j = i . The quantity 
k • is the volume average thermal conductivity between 

Wj 3 “i/i 

nodes j-1 and j, and k . . is the volume average thermal 

J + vl 

conductivity between nodes j and j+1. Figure 3a through d 
show the averages used for different interface locations. 

When X. = i k w< and k^ j+(/ = * s , which are the 

correct values to be used when the interface lies exactly 
on node j . It should also be noted that due to the method 
used to average k , the computer algorithm is, in general, 
only valid when the solid region is above the liquid region. 

Equations (41) through (43) are not used at the abrasion 
shield-ice interface or at the outer surface of the layer. At 
these surfaces, the boundary conditions (16) and (17) are 
finite differenced. The differencing procedure is essenti- 
ally the same as that used for the boundary conditions 
previously encountered. The finite difference equations 
used at these surfaces are given below: 
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1. Abrasion Shield - Ice Interface 

The grid shown in Figure 2a with i=II and i+b=w is used 
at this interface. Due to the dependence of k w on position, 
an alternate finite difference analog for k^ is used in 

equation ( 16b) . The analog is obtained by taking the average 
of equations (46a, b) ; that is: 

a . a , A 


‘wssr 


= — (k-» ^ Tw * k vr C ^ Tw l ^ (49) 


Then, substituting the above and analog (22) into boundary 
condition (16b) and finite differencing (16a) yields: 


i rp m 0 

IX, j - A w,j ' *11, j 


m O 

w,j 


(50a) 


*11 lEIdlL 


T 4- T 0 

H,j-1 H,j+1 


m O 

II,j-l 


4AX 


II 


k w, j+i/j T w, j+1 ~ T w, j* * k w, j-vi T w, j ~ T w, j-1* + 


(50b) 


4Ax w 

k° . (T ° . - T° ) 4- k° (T° - TO \ 

K w,j+|/ X u w,j+1 w,j ; XH 1 

Equation (26) with i=II is used to eliminate the fictitious 
temperatures T^ and j +1 , and equation (40) is used 
to eliminate the fictitious quantities k.„ .. (T 


T . . . ) and k° . (T ° . - T 0 . , 
w, j-1 w, j w, j-1 


w, D-l/x w, j 
) in (50b) . Combination 


of the resulting equation with equation (50a) and (15) 
yields. 



33 


for the solid range: 


Vj = ( H »°,j + Vw lT n,j-i + *A,j-i + < M n - »' 
T W °.j J + + N?,WV,j+l ” V.j )^ 1 + 

<M v/Ps e P s><Vj^i + “w + “iiV 1 


H . < H 
yf,3 smp 


(51a) 


with 


T . = H ./p£ 
w, 3 w , 3 f ns ps 


(51b) 


for the melting range: 


«w,j = H w,j + M wV 5 II,H - < M II + T mp + 


T II, j-1 + (M II " 1) T w° } j 1 + M w tk w, j+i/i T w, j+1 


- T ) + 
mp 


(52a) 


VO (T ° —TO \ 1 

W, j+‘/l W, j+1 w, j 7 1 


H < H . < H, 
smp w, 3 Imp 


with 


T . = T 
w ,3 mp 


(52b) 


for the liquid range : 


H w, j = (» w °, j + Vwl T n,H - <"ii + - 

+ *11, j-1 + (M XI “ l >V,jJ + Wj+^w.j+l 
+ “inp^pl 1 + k w. 3+i/i T w. j+1 - V.j'M 1 + 


(53a) 


- T 
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<VPiV <*„, j+1/1 + “w + *w 1 


(53a) 

continued 


H W,j S H 


Imp 


with 


T. 


w,D 


= (H, 


w, D 


- H 


)/pA 


lmp"HL pi 


+ T 


mp 


where 


\ = k ii 1 V 1K ii 


(53b) 


When node j is either solid or liquid, equation (51) or 
(53) reduces to equation (29). 

2 • Ice-Ambient Interface 

The grid shown in Figure 2d with II=w applies at this 
interface. Substituting the finite difference analogs (25) 
and (49) into boundary condition (17) gives: 


- k w, j+i/i T w, j+1 ~ T w, j) + k w, j-i/i T w, j ~ T w, j-1?* + 

4a *W 


V O (fp O . mO \ 4. V Q fT ° — T 0 ) 

Vj+i/lVj+i T w,r Vj-i/IVj 


(54) 


a2 


( i(T . + T 
\ 2 w, 3 


° . ) - T 
w,y a2 


) 


The quantities j+to <* w , j+1. ' T w,j> and k w, j+w (T w°, j+1 

- T° .) are eliminated using equation (40), and combination 
w , j 

of the result with (15) yields, for the solid range: 

H , = (h° . + M N 0 [2 T » - T° .] + M • 
w,D V w ,3 w a2,w a2 w , 3 w 




(55a) 
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( 1 ♦ <V'-°« 6 'p. ,<k w,j-'/ l + “a2,w>> 


(55a) 

continued 


with 


H w,j ^ H smp 

T w,j = H w,j^s C ps (55b) 


for the melting ranges 


V j = H w°, j + \ N a2,w C 2 T a2 - T mp " T w°,j J * 


(56a) 


H smp < H w, j < H lmp 


with 


T , = T 
v;, j mp 


(56b) 


for the liquid range: 

«w,l = ( H w 0 ,j + % “a2,w f 2 T a2 - V + H Xmp/ 


(57a) 


Pi e pl - T w°,j> + «w K,)-, ~ T mp + “imp/Pl^l* + 
k w,J-^ T w,j-l - V,j>^ 1 + <V-°L e pl )(k w,j-w + N a2,w>i 


H w, j > H lmp 


with 


where 


T w.j = <Vj - H lmp >/ -°i e pl + T , 


mp 


(57b) 


N a2,w ® h a 


Just like the other solid and liquid equations, (55) and 
(57) can be reduced to give equation (39) . 
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XV. COMPUTER IMPLEMENTATION 
A. METHOD OF SOLUTION 

The Crank -Nicolson finite differencing procedure 
results in atridiagonal set of linear equations that must 
be solved at each time level. Each equation relates the 
unknown temperatures T. . . T. . and T. ... at the current 
time level. The system of equations is solved by itera- 
tion because the phase of a node in the ice-water layer 
must be determined as the calculations proceed. The uae 
of matrix inversion methods would require the phase of a 
node to be known prior to the beginning of the calcula- 
tions. The Gauss-Seidel method was chosen as the iteration 
procedure because of its desirable convergence properties. 

The Gauss-Seidel method requires initial estimates of 
the value of the dependent variable at each node. These 
are obtained by either assigning the values calculated at 
the previous time step or by using linear extrapolation 
from the past to the present time level. Linear extrapo- 
lated values tend to speed up the convergence of the 
iteration. A series of passes is made through the grid, 
j=l,2,..., in which the value of the dependent variable is 
calculated at each node. This process is terminated when 
some convergence criterion is met. In these sweeps through 
the grid, the most current values are always used in the 
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calculations. For example, equation (26b) can be rewritten 
as: 


T i. ( r w = ( T i, < °“ > + T i, ( n > + T i°j+i + 2 <«i - 

T i,j + + 2 s i)/ 2 < M i + » 


(58) 


where is a value from the previous iteration and 

T . is a new value computed prior to computing t, 


All of the other finite difference equations can be 
rewritten in the same form as (58) . The convergence 
criterion used is that the difference between two succes- 
sive pass values must be less than some specified small 
value at each node. In most cases, 0.005$ was used. 

To accelerate convergence, which typically was slow, 
over-relaxation methods were used. These methods could 
not be applied to the ice-water layer equations due to 
stability* problems, and were only used for the composite 
blade. The successive over-relaxation (SOR) method yields 
the following modification of equation (53) : 


where T 


T (new) 

i, j 

(new) (58) 


i,3 


= (1 -U) T. ( ? ld > (59) 

*•>3 1,3 

is calculated from (58) . The parameter 


U) is known as the over-relaxation parameter, and acceler- 
ates convergence when 1<CJ < 2. The optimum value for U) 
varied from time level to time level, and was determined 
empirically. For the standard composite blade construction. 


it was found that the optimum CJ was about 1.7 for times 
less than 5 sec., 1.5 for times between 5 and 15 sec. and 

I. 3 for times greater than 15 sec. If the phase change is 
not considered, the ice layer can be treated as the II th 
layer of the blade, and over-relaxation is used for this 
layer also. 

The total number of calculations made can also be 
reduced by increasing the time step At as the calculations 
proceed. When ice shedding occurs, however, the rapid 
change in temperature Requires the standard time step of 
0.10 sec. to be reduced to 0.001 sec. in order for accurate 
results to be obtained. 

For more information on the methods used in the formu- 
lation and implementation sections, see References 10 and 

II . 

B . COMPUTER PROGRAM AND ALGORITHM 

The complete program listing appears in th« Appendix 
along with a sample input data file. The first eighty 
lines define the program input variables and their English 
units. The program can accept data in any consistent set 
of units, as only the input-output formats need be modified. 
A metric version of the program has been compiled and is 
available upon request. 

The flow chart for the main program is shown in Fig- 
ure 4. Dashed boxes may be skipped depending on the problem 
being solved. The subprograms used in the computer program 
have the following function: STEP determines the new time 
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step and adjusts time dependent parameters; SOURCE deter- 
mines the value of the source term at a half-time step; 
CONVE determines the percent difference between a new and 
old nodal temperature during an iteration; WLAYER calcu- 
lates the temperatures and enthalpies in the ice-water 
layer using the Enthalpy method; and PHASE determines the 

M. 

phase of a node and sets phase dependent properties. 

Figure 5 is the flow diagram for the subprograms WLAYER 
and PHASE, and illustrates the details of the determination 
of whether or not a given node is solid, liquid or melting. 
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V. DISCUSSION OP RESULTS 

The computer model developed in this investigation 
was used to study the effects of a number of design 
parameters on de-icer performance. These included the 
effects of heater power density and thickness, imperfect 
contact between layers, initial ice layer thickness, 
variable heater output, phase change in the ice layer and 
shedding of the ice. As stated previously, de-icer per- 
formance is measured by the time required to melt some 
specified thickness of ice at the abrasion shield-ice 
interface, starting from various initial composite blade 
temperatures. This time is referred to as the de-icing 
time; when it is reached, the ice can be shed by the 
dynamic forces acting on the outer surface of the composite 
blade. If the specified thickness of ice is zero, then the 
de-icing time is equal to the time required to raise the 
abrasion shield-ice interface temperature to 32®F. To deter- 
mine the zero thickness de-icing time, the phase change in 
the ice layer need not be considered. The phase change is 
not considered in Parts A through F below. Parts G through 
K require use of the Enthalpy method for the phase change. 

All figures and tables referred to below appear on pages 59 
through 85 . In addition to performance curves, temperature 
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response curves and profiles for some of the cases studied 
are presented and a comparison with experimental data is 
made. 

A schematic drawing of a composite aircraft blade on 
which an ice layer has formed appears in Figure la. Figure 
lb shows the one-dimensional finite difference grid used in 
the simulation and lists the number of nodes used in each 
layer for the standard de-icer design studied. These node 
numbers were enough to ensure that accurate solutions were 
obtained. Material property data and design data for the 
standard de-icer appear in Tables la and 2a, respectively. 
Any irariations from this design are clearly marked on all 
graphs presented. 


A. Verification of Finite Difference Method 

In order to verify the use of the Crank-Nicolson 
finite difference equations in the computer simulation, a 
problem for which an analytical solution existed was run. 
The problem chosen was to determine the temperature distri- 
bution in an infinite slab of thickness 2b as a function of 
time, when the slab was initially at a temperature T 0 and 
the surfaces of the slab were suddenly raised to a constant 
temperature T^. The analytical solution for the problem is 
[ 121 : 





( n+l/2 ) 


exp [ - ( n+l/2 ) 2 7 T 2 Qft/b 2 ] * 


cos [( n+1/2 my/bi 
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Th$ comparison between the analytical and numerical solutions 
appears in Figure 6. In this graph, the dimensionless tem- 
perature, T - T 0 is plotted as a function of the dimension- 



less distance, y/b, and the dimensionless time, Qft/b z . The 
analytical solution was obtained by summing a large number of 
terms in the infinite series. The comparison is very good, 
with only a slight discrepancy occurring for large times. A 
variable time step, which increased with time, was used in 
the simulation to achieve this accuracy. 

Baliga [1] also used the Crank -tlicolson finite difference 
scheme in his simulation. He compared the analytical and 
numerical solutions for a two layer slab problem and obtained 
equally good results. 

B. Effect of Power Density 

Figure 7 shows the effect of heater power density on 
de-icer performance. As with all performance graphs to 
follow, the temperature rise, which is the difference 
between the melting point of ice and the initial temperature, 
is plotted on the ordinate and the de-icing time is plotted 
on the abscissa. The power density curves computed ranged 
from 15 to 40 Watts/in 2 . These curves were also calculated 
by Baliga [1] and Stallabrass [2] with their computer models. 
There is perfect agreement between the results from this 
study and that of Baliga. This is expected since the phase 
change in the ice layer is not considered. Stallabrass' 
results tend to be slightly optimistic in comparison. The 
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curves shown in Figure 7 verify the observation made from 
experiments that the acceptable minimum power density is 
about 25 Watts/in 2 . The de-icing time increases rapidly 
as the power density is reduced, especially at low initial 
temperatures. 

C. Effect of Heater Thickness 

Heater thicknesses vary drastically depending on 
whether the heater is a woven mat of wires and glass fibers 
or a resistance ribbon. Woven mats are an order of magni- 
tude thicker than ribbons. Figure 8, however, shows that 
the heater thickness does not greatly affect de-icer per- 
formance, as the maximum difference between the de-icing 
times for the thicknesses shown is less than 1 sec. Curve 1 
is for a point (zero thickness) heater. This is an ideali- 
zation which shows the best possible results attainable. 
Curves 2 and 3 are for thicknesses characteristic of resis- 
tance ribbons, and curve 4 is for thicknesses characteristic 
of woven mats. 

D. Effect of Imperfect Contact between Layers 

The layers that make up a composite aircraft blade are 
held together by thin layers of epoxy resin. In addition, 
small air gaps may exist that cause poor contact between 
adjacent layers. These factors give rise to a resistance 
tc heat transfer across the layer interfaces. This resis- 
tance can be accounted for by means of an interfacial heat 
transfer coefficient. Figure 9 shows the effect of imper- 
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feet contact between layers on de-icer performance. The 
same heat transfer coefficient was assumed for all inter- 
faces in the blade, and perfect contact was assumed at the 
abrasion shield-ice interface. An infinite heat transfer 
coefficient corresponds to perfect contact. The results 
indicate that imperfect contact has little effect on per- 
formance down to a coefficient of about 500 Btu/ft 2 -hr-°F. 
Then, a drastic decrease in performance occurs between 500 
and 100 Btu/ft 2 -hr-°F, with de-icing times increasing by as 
much as 3 sec. The former coefficient corresponds to a 
resin thickness of approximately 0.005", and the latter to 
a thickness of approximately 0.01". This rather drastic 
change has not been accounted for in the previous investiga- 
tions surveyed . 

E. Effect of Initial Ice Layer Thickness 

In Figure 10, the effect of the ice thickness present 
when the heater is turned on is shown. The thicknesses 
studied ranged from 0.1" to 0.5". One might expect that 
the de-icing time would increase as the thickness of the 
ice layer is increased, but the opposite is true. The ice 
acts as a layer of insulation, so that the abrasion shield- 
ice interface temperature rises faster for thicker ice 
layers. For thin ice layers and high convection at the 
ice-ambient interface, the heat is rapidly conducted away 
from the abrasion shield-ice interface. In fact, the 
initial ice layer thickness has a greater effect on de-icer 
performance than any of the parameters previously discussed. 
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The curve for 0.1" of ice clearly shows that for low 
initial temperatures, thin layers of ice cannot be effec- 
tively removed with a power density of 25 Watts/in 2 . That 
the effect of the initial ice layer thickness has an upper 
limit can be seen by comparing the curves for 0.5" and 

O. 25" of ice. They are nearly identical. 

P. Effect of Variable Heater Output 

A variety of time-dependent heater outputs can be 
specified in the computer simulation. These are given in 
the "Numerical Formulation" section. A comparison was made 
between a heater with a sinusoidal output and a heater with 
a constant output of 25 Watts/in 2 , For simplicity, the 
phase change was not considered. The sinusoidal heater out- 
put that was studied is given by: 

q( t) =25 [ 1 + COS {IT/2 1 -TO 1 

It has an average power density of 25 Watts/in 2 and a period 
of 4 sec. The de-icing time starting from an initial tem- 
perature of -4°F was found to be slightly longer for the 
sinusoidal heater output, being 5.7 sec. compared to 4.9 
sec. for the constant heater output. 

The temperature responses for the substrate, heater, 
abrasion shield-ice interface and ice-ambient interface are 
shown in Figure ll for both heater outputs. For these 
response curves and those that follow in this section, the 
temperature variations across the substrate, heater and 
abrasion shield were usually much less than 1°F. The 
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sinusoidal heater output responses for the heater and the 
abrasion shield-ice interface oscillate, with different 
amplitudes and time lags, around the constant heater out- 
put responses. The substrate response, however, oscillates 
only slightly and is essentially super imposable with the 
constant heater output response. The temperature at the 
outer surface of the ice remains constant at -4°F. This 
is due to the large heat transfer coefficient at this 
surface . 

G . Application of the Enthalpy Method 

When the Enthalpy method is applied as described in the 
"Numerical Formulation" section, temperature responses like 
those shown in Figure 12 result. Figure 12 shows the abra- 
sion shield-ice interface temperature response for the 
standard de-icer design with a heater output of 25 Watts/in 2 
and an initial temperature of -4°F. The response behaves 
unrealistically after melting begins. Above 32° F, the tem- 
perature oscillates with a frequency which is nodal depen- 
dent. This can be seen by comparing the two curves for 20 
and 60 nodes in the ice layer. The broken curve is the 
temperature response predicted with the computer model of 
Baliga [1], and it also shows this behavior. These oscil- 
lations have been attributed to the fact that, when the 
Enthalpy method is used, a node in the ice layer remains 

at the melting point for a finite period of time. 

* 

Figure 12 shows that Baliga* s curve compares well in 
magnitude with those from the present study. The 3®F 
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difference which occurs after melting has begun is due to 
the fact that Baliga fl] approximates the latent heat effect 
with a large change in heat capacity over a 1° interval 
around 32 °P. The Enthalpy method does not require this 
approximation . 

Voller and Cross [8,9] have shown, by comparing ana- 
lytical and numerical solutions to simple phase change 
problems, that numerical solutions based on the Enthalpy 
method oscillate around the true solutions. They have also 
derived a criteria for determing the points of correspon- 
dence between the true and oscillating solutions. By finding 
these points of correspondence, accurate response curves can 
be obtained. It was shown in Section III, Part D, that the 
nodal enthalpy, H w ^ , could be directly related to the frac- 
tion of the node melted, Yy when the node was at the melting 

point. When Y.= l^ f the liquid-solid interface is exactly at 
3 2 

node j . It is when this occurs that the true and oscillating 

curves agree . By plotting the response variable when Y . = £ 

3 2 

at successive nodes, accurate response curves can be obtained. 

The above procedure was used to replot the 20 and 60 node 
curves in Figure 12, as well as to plot data for 30, 40 and 
90 nodes in the ice layer. The result is shown in Figure 13. 
The temperature response curve is now physically realistic 
and has very little nodal dependence. Thirty nodes was found 
to be the practical minimum number of nodes, and 30 nodes per 
every 0.25" of ice were used in all the results that follow. 
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Also, in all the graphs that consider melting in the ice 
layer, the points of correspondence obtained from the 
plotting procedure of Voller and Cross are clearly marked. 
The discontinuity in the slope which occurs when the abra- 
sion shield-ice interface reaches the melting point of 
32 ®F is characteristic of phase change problems. 

Also plotted in Figure 13 is the abrasion shield-ice 
interface temperature obtained with the approximation of 
equal liquid water and solid ice thermal conductivities. It 
is clear from Figure 13 that this is a bad assumption. In 
reality, liquid water has a lower thermal conductivity than 
ice. Thus, the thin layer of water which forms when the 
ice melts acts as an additional layer of insulation. This 
is the reason the true response curve lies above the approx- 
imate curve in Figure 13. 

It was found that the method of averaging used for the 

thermal conductivity between adjacent nodes in the ice layer 

significantly affected the numerical results. Only when a 

continuous volume average between adjacent nodes was used 

did the plotting procedure of Voller and Cross eliminate 

all nodal dependence. Volume averaging gives the correct 

values for k . and k . , , when the liquid-solid inter- 
w, 3-i/i w, 

face is exactly at node j (see Section III, Part D) . 

H. Effect of Phase Change 

The graphs discussed here are for the same set 
of conditions as were used for Figures 12 and 13. 
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Figure 14 shows the temperature responses for the 
substrate, heater, abrasion shield-ice interface and 
ice-ambient interface. Also plotted are the responses 
obtained without considering the phase change in the ice 
layer. The absorption of the latent heat during melting 
has the effect of lowering the temperatures in the composite 
aircraft blade after melting begins at 4.9 sec. 

The ratio of computing time to simulation time was 0.9 
without the phase change and 1.5 with the phase change on 
the University of Toledo IBM 4341 computer. For compari- 
son, Baliga's [1] computing times for the same problem 
were approximately five times longer (however, 60 nodes was 
the minimum number of nodes used in the ice layer) . This 
difference is partially due to the complexities of the 
matrix inversion technique used by Baliga. 

Figure 15 contains two temperature profiles across the 
composite blade-ice body. The profile for 4.0 sec. i3 
before melting begins, and the profile for 16.3 sec. is 
after a thin layer of ice has melted. The slight gradients 
across the substrate, heater and abrasion shield that were 
mentioned earlier are apparent. The profile after melting 
has begun contains an extra segment (#6) corresponding to 
the thin water layer which has formed next to the abrasion 
shield. 

The movement of the liquid-solid interface is plotted 
in Figure 16. The plotting procedure of Voller and Cross 
was also used for this response. For comparative purposes. 
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the curve for equal thermal conductivities and the curve 
predicted by Baliga's simulation [1] are shown. Baliga 
predicts a slightly longer abrasion shield-ice interface 
melting time than the present investigation, 5.5 sec. as 
compared to 4.9 sec. This is due to the approximate phase 
change technique used by Baliga. 

To directly consider the effects of heater element 
gaps on de-icer performance requires a two-dimensional 
model. This is because a finite thickness of ice will 
melt above the heater before any ice melts above the 
heater gaps. However, this effect can be studied indirectly 
by defining the de-icing time as the time required to melt 
different thicknesses of ice at the abrasion shield-ice 
interface.. This is done in Figure 17. The curve for 0" 
of ice is the same as the 25 Watts/in 2 <gnsrve in Figure 7. 

The other curves are for the different thicknesses of ice 
that must be melted for de-icing. They show the general 
increase in de-icing time with required thickness for 
de-icing. 

I. Comparison with Experimental Data 

Since the de-icer problems considered in this study 
cannot be solved analytically when the phase change is 
taken into account, a comparison of the numerical results 
obtained from the computer simulation with experimental 
data was made. Gent and Cansdale [3] present temperature 
data measured from three laboratory de-icer pads. The 
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material property data for these test specimens appears 
in Table lb. The design data appears in Tables 2b, c 
and d. 

Table 3 contains the numerical and experimental 
melting times for the abrasion 3hield-ice interface. 

The experiments were run for the three different de-icer 
pads at three different heater power densities, 16.6, 

19.0 and 22.5 Watts/in 2 . The initial temperature in all 
cases was 4.1°P. The ice layer was 0.1" thick and twelve 
nodes were used in this layer. It can be seen in column 3 
that When perfect contact between layers is assumed in the 
computer model, the predicted melting times are about 1 
sec. too short. Gent and Cansdale [3] give the thick- 
nesses for the glue between layers as between 0.001" and 
0.002". This enabled the interfacial heat transfer coef- 
ficient to be estimated at between 600 and 1200 Btu/ft 2 - 
hr-°P. These values were determined by dividing the thermal 
conductivity of epoxy resin by the glue thicknesses. The 
ranges for the melting times obtained using these heat 
transfer coefficients had a span of a few tenths of a 
second and appear in the fourth column of Table. 3. Almost 
all of the experimental data lies within these ranges . 

To further compare the numerical results from the 
computer simulation with experimental data, the abrasion 
shield temperature responses for all three specimens with 
16.6 Watts/in 2 were plotted. These curves are shown in 
Figure 18. Prior to the onset of melting the discrepancies 
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between the numerically predicted and experimental tem- 
perature data are less than 2°F. After melting begins, 
the temperature rises predicted by the computer model are 
too optimistic. The best comparison is for specimen 3. 
Considering the experimental difficulties encountered by 
Gent and Cansdale (Appendix of Reference 3), the agreement 
is still very good. 

J. Effect of Ice Shedding 

In actual operation, when the de-icing time is reached, 
the ice layer plus any water which has formed is shed by the 
dynamic forces acting on the outer surface of the blade. 
Then, once the heater has been turned off, a new ice layer 
may form. For this reason, the heater is turned on and off 
periodically. This process was simulated with the computer 
model by shedding the ice layer at the de-icing time and 
then, after a period of time, adding a new ice layer. The 
de-icing time was taken to be the time required to raise 
the abrasion shield-ice interface temperature to 32 °F. The 
ice was replaced every 20 sec., and the heater was turned 
cn for 10 sec. and then off for 10 sec. The temperature 
responses for the various locations in the composite air- 
craft blade are shown for the first 20 sec. cycle in 
Figure 19. A sharp decrease in the heater temperature 
occurs at 4.9 sec. and 10 sec. The first decrease is from 
the ice being shed, and the second is from the heater 
being turned off. The abrasion shield outer surface tem- 
perature drops immediately from 32® to -4 °F when the ice 
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is shed. Within 4 sec. of the heater being turned off , 
all the layers above the inner insulation have cooled to 
within 4°F of the initial temperature. Only the substrate 
remains hot after 20 sec. The time step had to be reduced 
drastically to follow the rapid change in temperatures 
which occurred when the ice was shed. 

The second cycle, from 20 to 40 sec., is shown in 
Figure 20. The responses are quite similar to those in 
Figure 19. The only significant difference is the sub- 
strate temperature which is hotter. The de-icing time 
decreased to 4.5 sec. The cyclic process was continued 
until a steady value of 4.4 sec. was obtained for the 
cyclic de-icing time. Thus, the difference between the 
first cycle de-icing time and the steady value is only 
0.5 sec. It is apparent that the temperature distribution 
present when the heater is turned on does not greatly 
affect de-icer performance. This is because once the ice 
is shed, heat is lost rapidly through the abrasion shield 
outer surface . 

K. Effect of Re freezing 

It is desirable to see what the temperature responses 
would be like if the ice layer could not be shed within the 
10 sec. heating time allotted in Part J. Figure 21 shows 
the temperature responses for this case. When the heater 
is turned off, the temperatures in the heater and abrasion 
shield drop immediately, with the heater temperature drop- 
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ping below the abrasion shield-ice interface temperature 
after 13.3 sec. Quite surprisingly, the water begins to 
refreeze at the abrasion shield outer surface after only 
14,5 sec. The computed temperature profile data reveals 
that this is not due to the complete refreezing of the 
water, but to the formation of a second ice layer. Thus, 
the water layer is sandwiched between two layers of ice . 
After 19.9 sec., the water has completely refrozen. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 

The one-dimensional electrothermal de-icer pad 
computer simulation developed in this study has been 
shown to predict accurate and consistent temperature 
distributions and ice-water interface location informa- 
tion. The accuracy was checked by comparing computed 
results with analytical solutions and experimental data. 
The simulation contains the following improvements over 
previous simulations on the subject: 

1) The Crank -Nicolson finite difference scheme was 
used instead of the forward finite difference 
scheme. This reduced the total number of calcu- 
lations that must be done by allowing a larger 
time step to be used. Round-off error was also 
decreased? 

2) The Enthalpy method was used for the phase change 
occurring in the ice layer. This method is more 
direct and easier to apply than methods previously 
used; 

3) The iterative method of solution used and the 
computer program algorithm can easily be extended 
to handle two-dimensional problems; and 

4) Many of the restrictions placed on previous simula- 
tions were removed. These restrictions included 
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perfect contact between layers In the composite 
aircraft blade and no consideration of the phase 
change in the ice layer. The present simulation 
can also consider a variety of boundary conditions, 
variable heater outputs and shedding of the ice 
layer. 

Further simulation work in the field of electrothermal 
de-icing should focus on the development of a two-dimen- 
sional computer model, so that the effect of heater gaps 
and blade geometry on de-icer performance can brs rigor- 
ously studied. This work is currently under way in the 
Chemical Engineering Department of the University of 
Toledo. In addition, a complete experimental study 
should be made on real de-icer pads. This would enable 
such parameters as heat transfer coefficients and layer 
thicknesses to be determined more accurately. It would 
also serve as a check on the assumptions used in the com- 
puter model. Finally, work should be initiated to 
determine and characterize the mechanism of ice deposi- 
tion on aircraft surfaces. 
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glue layers (10) 0.001 to 0.025 to between 1-2, 3-4, 4-5, 5-6 

0.002 0.05 
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Table 3. Comparison of Experimental and 
Predicted Melting Times 


Heater Time to Melt (see) 

Output 

(Watts/in 2 ) Simulation* Simulation* Experimental 


Specimen 1 

16.6 

7.9 

8. 6-8. 9 

9.0 


19.0 

6.7 

7. 3-7. 6 

7.5 


22.5 

5.4 

6. 0-6. 3 

6.3 

Specimen 2 

16.6 

5.8 

6.5-7. 0 

7.2 


19.0 

5.0 

5. 7-6. 2 

5.8 


22.5 

4.2 

4. 8-5. 2 

5.0 

Specimen 3 

16.6 

5.8 

6. 6-7.1 

7.3 


19.0 

5.1 

5. 7-6. 2 

5.9 


22.5 

4.2 

4. 8-5. 3 

5.0 


♦Assuming perfect contact between layers. 

+With contact resistance between layers. The lower value is 

Btu (0.001" of glue), and the higher 


for 


h. - 1200 

A* 


hr-ft 2 -°F 


value is for 


h L - 600 


Btu 


hr-ft 2 -°F 


(0.002" of glue). 


G 3P 
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Figure 2. Interface Finite Difference Grids 
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Figure 4. Flow Chart for Main Program 
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Figure 5. Flow Chart for Subprograms WLAYER and PHASE 
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De-Icing Time (sec.) 

Figure 7. Effect of Power Density on De-Icer Performance 


Standard De-Icer Design, 25 Watts/in 



De-Icing Time (sec.) 

Figure 8. Comparison of Point and Finite Thickness Heaters 



Figure 9. Effect of Imperfect Contact be 






Standard De-Icer Design, 25 Watts/in 



Figure 10. Effect of Initial Ice Layer Thickness on De-Icer Performance 







Standard De-Icer Design, 25 Watts/in 
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Figure 12. Abrasion Shield-Ice Interface Temperature with Phase Change 


Standard De-Icer Design, 25 Watts/in 





Figure 15. Temperature Profiles for Standard De-Icer Design Studied 


Standard De-Icer Design, 25 Watts 



Figure 16. Ice-Water Interface Location 






Figure 17. Effect of Ice Layer Thickness Melted on De-Icer Performance 



Experimental Test Specimens 1, 2 and 3, 16.6 Watts/in 



Figure 18. Comparison of Experimental and Predicted Abrasion Shield Temperatures 
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Figure 19. Temperature Response when Ice Layer is Shed, 1st Cycle 
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Figure 20. Temperature Response when Ice Layer is Shed., 2nd Cylce 


leer 



Figure 21. Temperature Response without Ice Shedding 
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Complete Program Listing 
and 

Sample Input Data Pile 


//IJuPT l 275 JOB ( UT , 

// 06200016, i) 

// EXEC PORT XC. LG 
//POUT . SYS IN DD * 

C 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


HEAT TRANSFER IN A COMPOSITE BODY 
INPUT VARIABLES: 

II, NUMBER OF LAYERS IN BODY, OR ONE LESS IF ICE-WATER LAYER IS 
INCLUDED. 

P1 = 0, NO HEATER; =1, VARIABLE WATTAGE HEATER; = 2, CONSTANT 
WATTAGE HEATER. 

P 2 = 0 , POINT HEATER; =1, FINITE THICKNESS HEATER. 

P3= I , FINITE THICKNESS HEATER; =2, POINT HEATER. 

P4 = 0, NONPERIODIC HEATER; =1, PERIODIC (ON-OFF) HEATER. 

P5=0, PHASE CHANGE NOT CONSIDERED; =1, PHASE CHANGE CONSIDERED. 
P6=0 , CONSTANT TEMPERATURE B.C. AT INNER SURFACE; =1, CONVECTIVE 
B.C. AT INNER SURFACE. 

P7 = 0, CONSTANT TEMPERATURE B.C. AT OUTER SURFACE; =1, CONVECTIVE 
B.C. AT OUTER SURFACE (P7.NE.0,IF PS=1). 

P8=0, CONSTANT TIME STEP USED; =1, VARIABLE TIME STEP USED. 

P9=0, NO PRINTING; =1, PRINT OUTPUT WHEN T(IX,JX) BECOMES .GE. 

TO TMAX (USED TO PRINT OUTPUT WHEN ICE BEGINS TO MELT) ; 

=2, TERMINATE PROGRAM AFTER ICE BEGINS TO MELT. 

P10=0, NO LINEAR EXTRAPOLATION; =1, LINEAR EXTRAPOLATION IS USED 
BETWEEN TIME STEPS TO ESTIMATE NEW TEMPERATURES. 

PI 1=0, CONSTANT ACCELERATION PARAMETER IS USED; =1, VARIABLE 

ACCELERATION PARAMETER IS USED FOR OVER-RELAXATION IN THE 
SINGLE PHASE LAYERS. 

P 1 2=0 , SHEDING OF ICE LAYER IS NOT CONSIDERED; =1, ICE LAYER 

IS SHED WHEN THE INTERFACE BETWEEN THE SOLID AND LIQ. IS 
AT NODE JSH ( ISM , DTMM AND DTMF SHOULD BE SPECIFIED). 

PI 3 = 0, INITIAL TEMPERATURE CONSTANT; =1, READ IN INITIAL 
TEMPERATURE DISTRIBUTION FOR COMPOSITE BODY. 

P14 = 0 , DO NOT STORE FINAL TEMPERATURES; =1, STORE FINAL 

TEMPERATURE DISTRIBUTION (P13 AND P.14 ARE USED FOR 
ICE SHEDING PROBLEMS) . 

JJ(I=1,II), NUMBER OF NODES IN LAYER. 

L (1-1,11) , THICKNESS OF LAYER (IN). 

K ( I = 1 , 1 1 ) , THERMAL CONDUCTIVITY OF LAYER ( BTU/HR-FT- ' F ) . 

DIF(I = 1,II), THERMAL DIFFUSIVITY OF LAYER (FT-FT/HR) . 

PI (1=1, II-l) =0, PERFECT CONTACT BETWEEN LAYERS I AND 1+1; =1, 

CONTACT RESISTANCE BETWEEN LAYERS I AND 1+1. 
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C H I ( I = 1 , 1 1-1 ) , HEAT TRANSFER COEFFICIENT BETWEEN LAYERS I AND 

C 1+1 (6TU/HR-FT-FT- 1 F) . 

C t Q r IF P2.EQ.0, HEATER IS LOCATED BETWEEN LAYERS IQ AND IQfl 

C (IQ.GE.l AND , LT .II); IF P2.EQ.1, HEATER IS LOCATED IN 

C LAYER IQ (IQ.GT.l AND . LT . I I ) . 

C C1,C2,C3,A4, A5, CONSTANTS IN THE EQUATION: 

C QF=C1*TM+C2+C3*C0S (A4*TM+A5) WHERE QF IS THE 

C TOTAL WATTAGE OF THE HEATER (WATTS/IN- IN) AND 

C TM IS TIME (SEC) . 

C TMO , LENGTH OF TIME THE (PERIODIC) HEATER IS ON (SEC). 

C TMF , LENGTH OF TIME THE (PERIODIC) HEATER IS OFF (SEC). 

C TIN , INITIAL TEMPERATURE OF BODY (*F). 

C TA1 , AMBIENT TEMPERATURE AT INNER SURFACE ( ' F) . 

C TA2 , AMBIENT TEMPERATURE AT OUTER SURFACE ( ' F) . 

C HI, HEAT TRANSFER COEFFICIENT AT INNER SURFACE ( BTU/HR-FT-FT- ' F) 

C H2, HEAT TRANSFER COEFFICIENT AT OUTER SURFACE (BTU/HR-FT-FT- ' F) 

C ISI, ISM, SEE BELOW. 

C DTMI , INITIAL TIME STEP (SEC). 

C DTMM , INTERMEDIATE TIME STEP, USED FOR TIME STEPS .GT.ISI AND 

C .LE.ISM (SEC). 

C DTMF , FINAL TIME STEP, USED FOR TIME STEPS .GT.ISM (SEC). 

C ISF, TOTAL NUMBER OF TIME STEPS. 

C IFRQ , FREQUENCY OF PRINTOUTS. 

C CMAX , CONVERGENCE CRITERIA FOR ITERATION OF TEMPERATURES (%) . 

C IWI , IWM , SEE BELOW. 

C WI, INITIAL ACCELERATION PARAMETER FOR OVER-RELAXATION. 

C WM, INTERMEDIATE ACCELERATION PARAMETER, USED FOR TIME STEPS 

C .GT.IWI AND .LE.IWM. 

C WF, FINAL ACCELERATION PARAMETER, USED FOR TIME STEPS .GT.IWM, 

C I X , JX , TMAX, SEE DESCRIPTION OF P9 . 

C CF1 , CONVERSION FACTOR FROM (IN-IN/SEC) TO (FT-FT/HR). 

C CF2 , CONVERSION FACTOR FROM (1/IN) TO (1/FT) . 

C CF3 , CONVERSION FACTOR FROM (WATTS/IN) TO ( BTU/HR-FT) . 

C CPS, SPECIFIC HEAT OF ICE (BTU/LB-'F). 

C KS, THERMAT CONDUCTIVITY OF ICE ( BTU/HR-FT- ' F) . 

C DENS, DENSITY OF ICE (LB/FT-FT-FT) . 

C MP, MELTING POINT OF WATER ( • F) . 

C DH, LATENT HEAT OF FUSION ( BTU/LB) . 

C CPL , SPECIFIC HEAT OF LIQ. WATER (BTU/LB-'F). 

C KL, THERMAL CONDUCTIVITY OF LIQ. WATER ( BTU/HR-FT- ' F) . 

C DENL, DENSITY OF LIQ. WATER (LB/FT-FT-FT) . 

C JW, NUMBER OF NODES IN ICE-WATER LAYER. 

C LW, THICKNESS OF ICE-WATER LAYER (IN). 

C JSH , SEE DESCRIPTION OF Pi 2. 

C 


o o o o o n 
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Aft? 


IMPLICIT REAL (K-N) 

INTEGER P1,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,PI,GS,TERM 
DIMENSION JJ(9) ,JN(9) ,L{9) ,DIF<9) ,K(9) ,DX(9) ,M(9) ,S(9) ,N{9) 
DIMENSION PI (9) ,HI(9) ,N1(9) ,N2(9) 

DIMENSION T { 70 ) ,TO(70) 

DIMENSION TW (91) ,TWO(91) ,H(91) ,HO(91) 

COMMON /AREA1/MP , HSMP , HLMP , HMP , CS , CL , TMP , KS , KL' 

COMMON /AREA2/JW , MW , NW , NA2W 
COMMON /AREA5/ISI, ISM,DTMM,DTMF 

COMMON /AREA6/P3,P4,CF3,TMO,PER,AA,Al,A2,A3, A4,A5 
DATA INS/56/, IN/5/, IO/6/, I0S/56/,CFl/25./,CF2/12./ 

CF3=40 .9463 

INPUT DATA FOR THE COMPOSITE BODY 

READ ( IN , 10) II,P1,P2,P3,P4,P5,P6,P7 
READ (IN,17)P8,P9,P10,P11,P12,P13,P14 
READ ( IN , 1 1 ) (JJ (I) ,1 = 1,11) 

READ ( IN , 12) ( L ( I ) , 1=1, II) 

READ ( IN , 1 2) ( K ( I ) ,1=1,11) 

READ (IN, 12) ( DIF ( I) , 1=1, II) 

III=II-1 

READ ( IN , 1 1 ) (PI (I) ,1=1,111) 

READ (IN, 12) (HI (I) ,1=1,111) 

READ (IN, 13) IQ,C1,C2,C3,A4,A5 
READ ( IN , 1 4 ) TMO , TMF 
READ (1N,15)TIN,TA1, TA2 , HI , H2 
READ (IN, 16) IS I , ISM,DTMI ,DTMM,DTMF 
READ (IN, 16) ISF,IFRQ,CMAX 
READ (IN, 16) IWI, IWM,WI,WM,WF 
READ ( IN , 16) IX, JX,TMAX 

INITIALIZATION AND CALCULATION OF TIME- INDEPENDENT CONSTANTS 

TM=0. 

IS = 0 
I P5=P5 
DTM=DTMI 
W=WI 
AA=0. 

PER=TMO+TMF 
TL=0 . 

TERM=0 
AGS=0 . 

GS = 0 
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JN(1) *JJ(1) 

DO 120 1=1,11 
TL=TL+L{ I) 

DX(I) -L(I)/(JJ (I) -1) 

M(I)=CF1*DX(I) **2/ (DIF (I) *DTM) 

S(I) =0. 

IF ( I . EQ, 1) GO TO 120 
JN(I)=JN(I-1)+JJ(I) 

IF (PI ( I — 1 ) . NE . 0 ) GO TO 116 

N(I-l) = K(I) *DX(I-1)/(K(I-1) *DX(I) ) 

GO TO 120 

116 N1 ( X- 1 ) =HI ( X-l ) *DX(I-1)/(K(I-1) *CF2) 

N2(I-1)=HI (1-1) *DX(I)/(K(I) *CF2) 

120 CONTINUE 

IF ( P9 . EQ. 0 ) GO TO 122 
JXX=JX 

IF ( IX. NE. 1) JXX=JXX+JN ( I X-l ) 

122 IF (P13.EQ.0)GO TO 124 

INPUT INITIAL TEMPERATURE DISTRIBUTION FOR COMPOSITE BODY 


JN0 = 1 

DO 123 1=1,200 
JNl=JN0+5 

IF ( JN1 .GT . JN (II)) JN 1=JN (II) 
READ ( INS, 18) (T ( J) ,J=JN0,JN1) 
IF (jNl . EQ. JN (II) ) GO TO 127 
JN0=JN1+1 

123 CONTINUE 

124 JN1= JN (II) 

DO 125 J ~1 , JNl 
V(J) =TIN 

125 CONTINUE 

127 IF ( P6 . EQ. 0) GO TO 128 
NA1=H 1 *DX ( 1 ) / ( K ( 1 ) *CF2) 

GO TO 130 

128 T ( 1 ) =TA1 

130 IF ( P7 - CQ . 0) GO TO 132 

NA2=112*OX ( 1 1 ) / (K ( 1 1 ) *CF2) 

GO TO 135 

132 IF(P5.EQ.0) T ( JN ( I I ) ) =TA2 
135 IF ( PI * EQ . 0 ) GO TO 140 
IF (P2 . EQ. 0 . ) GO TO 138 
A1=C 1 *DX (IQ) /L { IQ) 

A2=C 2*DX (IQ) /L( IQ) 
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138 


A3=C3*DX(IQ)/L(IG>) 
GO TO 140 


Al =01 
A2=C2 
A3-C3 


op'S «« IS 
0F P0 °R quality 


PRINT DATA FOR THE COMPOSITE BODY 


140 WRITE (10,20) 

WRITE (10,21) 

WRITE (10,22) 

IF ( P5 . EQ . 0 ) GO TO 145 
IPRT-II+1 
WRITE (10,23) I PRT 
WRITE (10,24) 

GO TO 150 
145 WRITE (10,23) II 
WRITE (10,25) 

150 WRITE (10,26) 

WRITE (10,27) 

WRITE (10,28) 

WRITE (10,29) 

WRITE (10,30) 

DO 160 1=1,11 

WRITE ( 10, 31) I, JJ (I) , L ( I ) , DX (I) , K ( I ) , DI F ( I ) 
160 CONTINUE 

WRITE (10,32) TL 
I F ( P 1 . EQ . 0 ) GO TO 200 
WRITE (10,33) 

WRITE (10,34) 

IF (P2 . EQ. 1 . ) GO TO 170 
IQ1 =IQ+1 

WRITE (10,35) IQ , IQ1 
GO TO 180 
170 WRITE (10,36) IQ 
180 IF ( P4 . EQ . 0 ) GO TO 190 
WRITE(IO,37) TMO 
WRITE (10,39) TMF 
190 WRITE (10,40) 

WRITE (10,41) 

WRITE(I0,42)C1,C2,C3 
WRITE(I0,43) A4,A5 
200 WRITE (10,44) 

WR ITE (10,45) 

WRITE ( 10 , 4 6 ) TIN 
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I F ( P6 , NE . 0 ) GO TO 210 
WRITE (10, 47) TA1 
GO TO 220 
210 WRITE (10,48) 

WRITE ( 10 , 49 } TAl 
WRITB(IO, 50) HI 
220 IF (P7 .NE. 0) GO TO 230 
WRITE(I0,51) TA2 
GO TO 240 
230 WRITE (10,52) 

WRITE ( 10 , 4 9 ) TA2 
WRITE ( 10, 50) H2 

240 DO 241 1=1,111 

IF (PI (I) . EQ . 0 ) GO TO 241 
IPRT=I+1 

WRITE (10, 87) I, IPRT 
WRITE (10, 50) HI (I) 

241 CONTINUE 
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INPUT DATA FOR ICE-WATER LAYER IF INCLUDED 


IF (P5 . EQ . 0) GO TO 260 
READ (IN, 14) CPS, KS, DENS, MP,DH 
READ ( IN , 14 ) CPL, KL, DENL 
READ ( IN , 16 ) JW , JSH , LW 

INITIALIZATION AND CALCULATION OF TIME-INDEPENDENT CONSTANTS 

HSMP=DENS*CPS*MP 
HLMP=DENL* (CPS*MP+DH) 

HMP= . 5* ( HLMP+HSMP) 

CS=DENS*CPS 

CL=DENL*CPL 

TMP=MP-HLMP/CL 

IF (TIN . LE.MP) HIN=CS*TIN 

IF(TIN.GT.MP) HIN=CL* (TIN-TMP) 

DO 247 J=1,JW 
TW ( J ) =TIN 
H(J) =HIN 
247 CONTINUE 

DXW=LW/ ( JW-1 ) 

MW=DTM/(DXW**2*CF1) 

NW=K (II) *DXW/DX (II) 

NA2W=H2*DXW/CF2 
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PRINT DATA FOR THE ICE-WATER LAYER 


WRITE (10, 53)) 

WRITE (10,54) 

WRITE (10,55) 

WRITE (10,56) DENS , DENL 
WRITE ( 10 , 57 ) CPS , CPL 
WRITE (10,58) KS, KL 
WRITE (10, 59) HSMP,HLMP 
WRITE (10, 60)MP 
WRITE (10,61) DH 
WRITE ( 10 , 62 ) TIN 
WRITE (10,63) HIN 
WRITE (10, 64)1LW 
WRITE(IO,65 ) JW , DXW 
IF (P12.EQ.0)GQ TO 260 
WRITE ( 10 , 96 ) JSH 

260 WRITE (10, 66) 

WRITE (10,67) 

WRITE (10,95) IN 

IF ( P8 . EQ. 0 ) GO TO 270 

WRITE (10,68) 

WRITE (10,69) DTMI , IS I 
ISI1=ISI+1 

WRITE (10,70) DTMM, IS II , ISM 
ISMl = ISM-fl 

WRITE (I 0 , 7 1 ) DTMF , ISMl 
GO TO 280 

270 WRITE(IO,72) DTMI 

280 WRITE (10,73) IFRQ 
WRITE (10,74) CM AX 
IF ( PI 0 . EQ . 0 ) GO TO 282 
WRITE (10,86) 

282 IF ( PI 1 . EQ . 0 ) GO TO 284 
WRITE (10,90) 

WRITE ( 10, 91) WI , IWI 
IWI1=IWI+1 

WRITE (10,92 )WM, IWI1, IWM 
IWM1= IWM+1 
WRITE ( 10 , 9 3 ) WF , IWM1 
GO TO 286 

284 WRITE ( 10,94 ) WI 

286 WRITE ( 10, 75) 

WRITE (10,76) 

WRITE( 10,77) 
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IF (P 1 3 . EQ . 0 ) GO TO 290 
WRITE (£0*78) TM,GS 
JN0 = 1 

DO 287 1*1,11 
WRITE (10,79) I 
JNi-JN { I ) 

WRITE (10,80) <T(J) ,J=JN0,JN1) 

JN0=JN1 .+ 1 
I. CONTINUE 

IF(P5.EQ.0)GO TO 290 

I PRT-I 1+1 

WRITE (10,81) I PRT 

WRITE (10,80) ( TW ( J ) ,J=1,JW) 

WRITE (10,82) 

WRITE (10/89) (H (J) ,J=1,JW) 

INITIALIZATION OF VARIABLES FOR NEW TIME STEP 

290 IS=IS+1 
P5=IP5 
JN1=JN (II) 

DO 300 J=1,JN1 

I F ! { PI 0 . EQ . 0 . OR . IS . EQ . 1 ) GO TO 292 

LINEAR EXTRAPOLATION 

TEXT«?2 . *T ( J) -TO ( J) 

292 TO ( J) =T ( J) 

IF(P10.EQ.0. OR . IS . EQ . 1 ) GO TO 300 
T ( J ) =TEXT 
300 CONTINUE 

IF(P5.EQ.0)GO TO 320 
DO 310 J=1,JW 

IF(P10.EQ.0. OR .IS.EQ.1)G0 TO 302 

LINEAR EXTRAPOLATION 

HEXT=2 . *H ( J ) -HO ( J) 

302 HO ( J ) =H ( J) 

TWO ( J ) =TW ( J) 

IF(P10.EQ.0.OR.IS.EQ.1)GO TO 310 
H (J) =HEXT 

I F ( H ( J ) .LE.HSMP) JP=1 

I F ( H ( J ) .GT.HSMP.AND.H (J) .LT.HLMP) JP=2 
I F ( H ( J ) .GE.HLMP) JP=3 
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GO TO 

(304,305, 

306) ,JP 

304 

TW ( J ) 

=H(J)/CS 



GO TO 

310 


305 

TW ( J ) 

=MP 



GO TO 

310 


306 

TW ( J ) 

=H ( J ) /CL4 

■TMP 

310 

CONTINUE 


320 

IF ( P8 

. EQ . 0 ) GO 

TO 330 


ADJUSTMENT OF CONSTANTS IF TIME STEP CHANGES 

CALL STEP (P5 , I I , IS , DTM,M,MW) 

330 IF(P11.EQ.0) GO TO 332 

ADJUSTMENT OF ACCELERATION PARAMETER 


IF ( IS . EQ . IWI+1 ) W=WM 
IF ( IS. EQ . IWM+1 ) W=WF 
332 IF(P1.EQ.0) GO TO 340 

I F ( P 1 . EQ . 2 . AN D . I S . CT « 1 ) GO TO 340 

CALCULATION OF NEW SOURCE TERM 

CALL SOURCE ( DX (IQ) ,K(IQ) ,S(IQ) ,TM,DTM) 

CALCULATION OF NEW TIME 

340 TM=TM+DTM 

GAUSS -SEIDEL REITERATION 

DO 400 GS=1 ,200 
ICV=0 

IF ( P6 . EQ. 0 ) GO TO 350 

CALCULATION OF TEMPERATURE AT THE INNER-AMBIENT INTERFACE 
TOLD=T ( 1 ) 

T ( 1 ) =T { 1 ) +W* ( (T(2)+TO(2)-( (X.-M(l) )+NAl) *TO(l) 

1+2 . *NA1*TA1 ) / ( 1 . +M ( 1 ) +NA1 ) -T ( 1 ) ) 

CALL CONVE(TOLD,T(l) ,CMAX, ICV) 

350 JN0=2 

DO 360 I=l r II 
JN 1= JN (I ) — 1 

IF ( JN1 . LT . JN0) GO TO 356 
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CALCULATION OF TEMPERATURES IN THE INTERIOR OF THE LAYER 
TOLD=T(J) 

T(J)=T(J)+W*( (T(J+l)+T(J-l)+TO(J+l)+2.*(M(I)-l.)*TO(J) 

1+TO ( J-l ) +2 . *P2*S ( I ) ) / ( 2 . * (M ( I ) +1 . ) ) -T(u) ) 

IF(ICV.riE.0)GO TO 355 

CALL CONVE (TOLD , T (J ) ,CMAX,ICV) 

355 CONTINUE 

356 IF ( I . EQ . I I ) GO TO 370 

IF ( PI ( I ) . NE . 0 ) GO TO 357 

CALCULATION OF TEMPERATURE AT THE INTERFACE BETWEEN LAYERS, 
NON-RESISTIVE INTERFACE 

TOLD=T (JNl+1 ) 

T(JNl + .l)=T (JNl + 1) -W* (T(JN1 + 1)-(T(JN;1)+N(I) *T ( JN1+3) +TO ( JN1) 
l+( (M(I)-1.)+N(I) *(M(I+1)-1.) ) *TO ( JNi + 1 ) +N ( I ) *TO ( JN1 + 3 ) 

1+S ( I ) +P2*N ( I ) *S(I+1) ) / ( 1 . +M ( I ) +N ( I ) * ( 1 * +M ( I + 1 ) ) ) ) 

T { JN1 + 2) =T (JNl+1) 

IF ( ICV . NE . 0 ) GO TO 359 

CALL CONVE (TOLD,T ( JNI+I) ,CMAX,ICV) 

GO TO 359 

CALCULATION OF TEMPERATURES AT THE INTERFACE BETWEEN LAYERS, 
RESISTIVE INTERFACE 

357 TOLD=T ( JNl+i ) 

T(J , U+X)“T(JN1+1)+W*( (T(JN1)+N1(I) * (T ( JN1+2) +TO ( JN1+2) ) - 
1 (l.-M(I)+Nl (I) ) *TO (JNl+1 ) +TO ( JN1 ) +S ( I ) ) / ( 1 . +M ( I ) +N1 ( I ) ) -T (JNl+1) ) 
IF ( ICV. NE . 0 ) GO TO 358 
CALL CONVE (TOLD, T (JNl+1) ,CMAX,ICV) 

358 TOLD=T ( JN1+2) 

T(JNl+2)=T(JNl+2) +W* ( (T(JNl+3)+N2(I) * (T (JNl+1 ) +TO (JNl+1 ) }- 
1(1. -M ( 1 + 1 ) +N2 ( I ) ) *TO ( JN1+2) +TO ( JN1 + 3) +P2*S ( 1 + 1 ) ) / 

1(1. +M { 1 + 1 ) +N2 ( I ) ) -T ( JN1+2) ) 

IF ( ICV . NE . 0 ) GO TO 359 

CALL CONVE (TOLD, T(JNl+2) ,CMAX,ICV) 

359 JN0=JNl+3 

360 CONTINUE 

370 IF (P5.NE, 0) GO TO 385 
IF ( P7 . EQ . 0 ) GO TO 390 

CALCULATION OF TEMPERATURE AT THE OUTER-AMBIENT INTERFACE 
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C 

TOLD=T ( JN1 + 1 ) 

T(JN1+1)=T(JN1+1)+W*( (T(JNl)+TO(JNl)-( (l.-M(II) )+NA2) 
1*T0(JN1+1)+2.*NA2*TA2)/(1.+M(II)+NA2) -T(JN1+1) ) 

IF ( ICV. NE. 0) GO TO 390 

CALL CONVE(TOLD,T(JNl+l) ,CMAX,ICV) 

GO TO 390 
C 

C CALCULATION OF TEMPERATURES IN THE ICE-WATER LAYER 
C IF INCLUDED 
C 

385 CALL WLAYER (TW , TWO, H , HO, T ( JN1 ) ,TO(JNl) ,T(JN1+1) ,TA2, 
1M(II) , CMAX , ICV , GS) 

C 

C CHECK CONVERGENCE OF THE GAUSS-SEIDEL ITERATION 
C 

390 I F ( ICV . EQ . 0 ) GO TO 410 
400 CONTINUE 
410 WRITE (10,78) TM,GS 
AGS=AGS+FLOAT (GS) 

C 

C DETERMINATION OF WHETHER ICE LAYER SHOULD BE SHED 
C 

IF(P5.EQ.0. OR ,P12.EQ.0)GO TO 415 
HJSH=H ( JSH) 

IF ( JSH . EQ , 1) HJSH=HMP+ . 5* ( H ( JSH) -HSMP) 

IF (HJSH . LT. HMP) GO TO 415 
IP5=0 
P8 = 1 
ISI=IS 
C 

C DETERMINATION OF WHETHER THE PROGRAM SHOULD BE TERMINATED 
C 

415 IF(IS.EQ.ISF) TERM=1 
C 

C DETERMINATION OF WHETHER OUTPUT SHOULD BE PRINTED 
C 

IF(P5.NE.IP5)GO TO 430 
I F ( P9 . EQ . 0 ) GO TO 420 

IF(TO(JXX) . LT . TMAX . AND . T ( JXX) .GE.TMAX)GO TO 429 
420 IF (IS/IFRQ*I FRQ . EQ . IS . OR . TERM. NE . 0 ) GO TO 430 
GO TO 290 
C 

C PRINT OUTPUT OF PROGRAM 

C 
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429 

430 


440 


455 


STORE TEMPERATURE DATA FOR NEXT RUN 

460 IF ( P14 . EQ . 0 ) GO TO 465 
JN0 = 1 

DO 462 1=1,200 
JN1-JN0+5 

IF ( JN1 .GT . JN (II) ) JN1=JN (II) 

WRITE ( IOS , 18 ) (T ( J ) ,J=JN0,JN1) 

IF(JN1.EQ.JN (II) )GO TO 465 
JN0=JN1+1 
462 CONTINUE 

PRINT REASON WHY PROGRAM WAS TERMINATED 

465 WRITE (10,85) IS 

AVER=-AGS/FL0AT ( IS ) 

WRITE(I0,88) AVER 
STOP 

FORMAT STATEMENTS FOR INPUT AND OUTPUT 

10 FORMAT (3X,I2,5X,I1,5X,F1.0,5X,F1.0,5X,I1,5X,I1,5X,I1,5X,I1) 

11 FORMAT (5X, 818) 

12 FORMAT (5X,8F8.0) 

13 FORMAT (5X,I6,5X,F6.0,5X,F6.0,5X,F6.0,5X,F6.0,5X,F6.0) 

14 FORMAT (5X,F6.0,5X,F6.0,5X,F6.0,5X,F6.0,5X,F6.0,5X,F6.0) 

15 FORMAT (5X,F6.0,5X,F6.0,5X,F6.0,5X,F8.0,5X,F8.0) 

16 FORMAT (5X,I6,5X,I6,5X,FC.0,5X,F6.0,5X,F6.0,5X,F6.0) 


IF ( P9 . EQ. 2) TERM=1 
JN0 = 1 

DO 440 1=1,11 
WRITE (10,79) I 
JN 1=JN ( I ) 

WRITE (10,80) (T ( J ) , J =JN0 , JN 1 ) 

JN0=JN1+1 

CONTINUE 

IF ( P5 . EQ. 0 ) GO TO 455 

I PRT=I 1+ 1 

WRITE (10,81) I PRT 

WRITE (10,80) (TW ( J ) , J=1 , JW) 

WRITE (10,82) 

WRITE (10,89) (H ( J ) ,J = 1,JW) 

IF (TERM . NE . 0 ) GO TO 460 
GO TO 290 
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17 FORMAT (10X,I1,5X,I1,5X,II,5X,I1,5X,I1,5X,I1,5X,I1) 

10 FORMAT ( 6F1 2 . 5 ) 

20 FORMAT ( 1 9X , ' ' ) 

21 FORMAT ( 19X , ' HEAT TRANSFER ANALYSIS OF A COMPOSITE BODY') 

22 FORMAT ( 1 9X , ' ') 

23 FORMAT (// , 5X , ' THERE ARE % 1 2 , X , ' LAYERS IN THE BODY.*) 

24 FORMAT (/ , 5X, ' THE PHASE CHANGE IN THE ICE-WATER«r£jAYER 
1 IS CONSIDERED. ' ) 

25 FORMAT */,5X, 'A PHASE CHANGE IS NOT CONSIDERED.') 

26 FORMAT (//, 5X, ' PHYSICAL PROPERTY DATA FOR LAYERS') 

27 FORMAT (5X, ' •) 

28 FORMAT(/,5X, 'LAYER' ,4X, 'NUMBER' ,5X, 'LENGTH' ,6X, 'DELTA X' ,6X, 

1 'THERMAL' ,7X, 'THERMAL' ) 

29 FORMAT(5X, 'NUMBER' ,2X, 'OF NODES ', 26X , ' CONDUCTIVITY 3X , 

1 ' DIFFUS I VITY ' ) 

30 FORMAT ( 26X , * (IN) ' ,8X, • (IN) ' ,4X, ' (BTU/HR-FT- 1 'F) ' ,2X, ' (FT-FT/HR) ' 
1) 

31 FORMAT */,6X,I2,8X,I2,6X,F7.4,6X,F7.5,6X,F6.3,8X,F6.4) 

32 FORMAT */,9X, 'TOTAL LENGTH = ' f F7.4) 

33 FORMAT ( // , 5X , ' DATA FOR HEATER') 

34 FORMAT (SX, ' ') 

35 FORMAT */,5X, 'THE HEATER IS A POINT HEAT SOURCE LOCATED 
1 BETWEEN LAYERS' ,I2,X, 'AND' ,12, ) 

36 FORMAT (/,5X, 'THE HEATER HAS A FINITE THICKNESS , AND IS 
1 LOCATED IN LAYER ',12,'.') 

37 FORMAT*/, 5X, 'THE HEATER IS ON PERIODICALLY 10X ,' TIME ON - ', 
1F6 . 3 , X, * (SEC) ' ) 

39 FORMAT*/, 45X, 'TIME OFF = ' ,F6. 3,X, ' (SEC) ’ ) 

40 FORMAT (/,5X, 'THE TOTAL HEAT GENERATION IS GIVEN BYt') 

41 FORMAT (/, 1 3X, *QF = C1*TM + C2 + C3*COS(C4*TM + C5) ', 5X ,' (WATTS/ 
1 IN-IN) ' ) 

42 FORMAT (/ , 1 3X , 'Cl = ' ,F7 . 3 , 10X, ' C2 = ' , F7 . 3 , 1 0X , ' C3 = ',F7.3) 

43 FORMAT (/ , 1 3X, ' C4 = ' , F7 . 3 , 10X, ' C5 = ',F7.3) 

44 FORMAT*//, 5X, ' INITIAL AND BOUNDARY CONDITIONS') 

45 FORMAT (5X, ' ') 

46 FORMAT*/, 5X, ' INITIAL TEMPERATURE IN THE COMPOSITE BODY = 

1F6 . 2 , X , ' ( ' ' F) ' ) 

47 FORMAT (/,5X, 'CONSTANT TEMPERATURE AT THE INNER-AMBIENT 
1 INTERFACE = ' , F6 . 2 , X , ' ( " F) ' ) 

48 FORMAT (/,5X, 'CONVECTIVE BOUNDARY CONDITIONS AT THE INNER- 
1 AMBIENT INTERFACE:') 

49 FORMAT*/, 19X, 'AMBIENT TEMPERATURE = ' , F6 . 2 , X , ' ( " F) ' ) 

50 FORMAT*/, 13X, 'HEAT TRANSFER COEFFICIENT = ',F10.2,X, 

1 ’ ( BTU/HR-FT-FT- ' ' F) • ) 

51 FORMAT */,5X, 'CONSTANT TEMPERATURE AT THE OUTER-AMBIENT 
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1 INTERFACE = ' , F6 , 2 , X , ' ( ’ ' F) ' ) 

52 FORMAT*/, 5X, ’CONVECTIVE BOUNDARY CONDITIONS AT THE OUTER- 
1AMBIENT INTERFACE:') 

53 FORMAT *//,5X, 'PROPERTY DATA FOR THE ICE-WATER LAYER') 

54 FORMAT (5X, ' — ' ) 

55 FORMAT*/, 5X, ’PROPERTIES OF WATER :', 1 3X ,' ICE ’, 8X ,' LIQ . 

I WATER') 

56 FORMAT */,5X, 'DENSITY' , 18X , ' = ' , 4X , F8 . 4 , 6X , F8 . 4 , 2X , 

1' (LBS/FT-FT-FT) ’ ) 

57 FORMAT*/, 5X, 'SPECIFIC HEAT' , 1 2X , ' - ' , 4X , F8 . 4 , 6X , F8 . 4 , 2X 
1' ( BTU/LB- ' ' F) ' ) 

58 FORMAT*/, 5X, 'THERMAL CONDUCTIVITY' , 5X , ' = ' , 4X , F8 . 4 , 6X , F8 . 4 , 2X , 

1' ( BTU/HR-FT- ' * F) ’ ) 

59 FORMAT*/, 5X, ' ENTHALPY AT MELTING PT . ' , 2X , ' - ' , 3X , F9 . 4 , 5X , F9 . 4 , 2X 
1’ * BTU/FT-FT-FT) ' ) 

60 FORMAT */,5X, 'MELTING POINT OF WATER ' , 10X , ' - ' , 8X , F5 . 2 , X , ' ( ' ' F) ' ) 

61 FORMAT*/, 5X, ' LATENT HEAT OF FUSION ', 11X ,'=', 5X, F8 . 3 , X ,' (BTU 
1/LB) ' ) 

62 FORMAT*/, 5X, ' INITIAL TEMPERATURE IN LAYER ' , 4X , ’ = ' , 7X , F6 . 2 , X , 

1* ( « VF) ' ) 

63 FORMAT*/, 5X, ' INITIAL ENTHALPY IN LAYER ' , 7X , ' = ' , 4X , F9 . 3 , X, 

1' (BTU/FT-FT-FT) ' ) 

64 FORMAT*/, 5X, ’LENGTH OF LAYER ' , 17X , ' = ’ , 6X, F7 . 4 , X, ' * IN) ! ) 

65 FORMAT */,5X, 'THERE ARE ' , 1 3 , X , ' NODES AND DELTA X ’ , 2X , ' - ' , 6X , 

1F7 . 5 , X , ' (IN) ' ) 

66 FORMAT*//, 5X, 'ADDITIONAL DATA’) 

67 FORMAT (5X, ' •) 

68 FORMAT*/, 5X, 'A VARIABLE TIME STEP IS USED:') 

69 FORMAT*/, 5X, ' INITIAL TIME STEP ', 6X, X, *■> >, X ,’ (SEC) ,' , 

1' TIME STEPS 1 THROUGH', 14) 

70 FORMAT*/, 5X, ' INTERMEDIATE TIME STEP = ' ? F4 . 2 , X , ' (SEC) , ' , 

1' TIME STEPS' ,14, X, 'THROUGH' ,14) 

71 FORMAT*/, 5X, 'FINAL TIME STEP ' , 8X , ' * ' , X, F4 . 2 , X , ' * SEC ) , ' , 

1* TIME STEPS' ,14, X, 'ON' ) 

72 FORMAT */,5X, 'CONSTANT TIME STEP = * , X , F4 . 2 , X , ' ( SEC) ' ) 

73 FORMAT*/, 5X, 'OUTPUT IS PRINTED EVERY' , 12, X, 'TIME STEPS.') 

74 FORMAT */,5X, 'THE CONVERGENCE CRITERIA FOR TEMPERATURE 
1 IS ' ,F6.4,X, '%. ' ) 

75 FORMAT*//, 24X, ' ' ) 

76 FORMAT(24X, 'TEMPERATURE PROFILES IN DEGREES F' ) 

77 FORMAT ( 24X , ' ') 

78 FORMAT*//, 24X, 'TIME =' ,F8 . 3,X, • (SEC) , ' ,5X, 'GS = ' ,13) 

79 FORMAT*/, 5X, ' LAYER' , 12) 

80 FORMAT*/, 6F12. 5) 

81 FORMAT*/, 5X, ' LAYER' , 12, ' : ICE-WATER LAYER') 
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82 FORMAT ( / , 5X , 'ENTHALPY IN ICE-WATER LAYER, BT’J/FT-FT-FT' ) 

85 FORMAT (// , 5X , 1 THE PROGRAM WAS TERMINATED AFTER ', 14 , X, ' TIME 
1 STEPS WERE COMPLETED.') 

86 FORMAT (/, 5X, ‘LINEAR EXTRAPOLATION IS USED BETWEEN TIME STEPS'// 
17X, 'TO INITIALIZE TEMPERATURES.') 

87 FORMAT (A 5X, 'INTERFACIAL RESISTANCE BETWEEN LAYERS' 

1 , 1 2 , X , ' AND ' ,12,' :') 

88 FORMAT (/, 5X, 'THE AVERAGE NUMBER OF ITERATIONS PER TIME 
1STEP WAS' , F6 . 2 , X , ' . ' ) 

89 FORMAT (/, 6F12. 3) 

90 FORMAT (/,5X, 'A VARIABLE ACCELERATION PARAMETER IS USED FOR 
1 OVER-RELAXATION: ' ) 

91 FORMAT(/,5X, * INITIAL PARAMETER ', 6X , ' = ' ,X, F4 . 2 , X, • TIME 
1 STEPS 1 THROUGH ',14) 

92 FORMAT(/,5X, ' INTERMEDIATE PARAMETER - ' , X , F4 . 2 , X , ' , ' , • TIME 
1 STEPS' ,14, X, 'THROUGH' ,14) 

93 FORMAT(/,5X, 'FINAL PARAMETER' ,8X, ' =» ,X,F4 . 2, X, • TIME 
1 STEPS' ,I4,X, 'ON' ) 

94 FORMAT (/,5X, 'THE CONSTANT ACCELERATION PARAMETER FOR 
1 OVER-RELAXATION IS ' , X , F4 . 2 , X , ' . ' ) 

95 FORMAT(/,5X, ' INPUT DATA FILE ',I3,X,'.') 

96 FORMAT (/, 5X, ' LAYER IS SHED WHEN THE ICE-WATER INTERFACE IS 
1 AT NODE' , 13, ' . ' ) 

END 

SUBROUTINE STEP ( P5 , I I , IS , DTM , M , MW) 

STEP DETERMINES NEW TIME STEP AND ADJUSTS TIME-STEP DEPENDENT 
CONSTANTS 

REAL M ( 6 ) , MW 

COMMON /AREA5/IS I , ISM , DTMM , DTMF 
A=DTM 

IF ( IS . EQ. ISI+1) DTM=DTMM 
IF ( IS . EQ . ISM'f 1 ) DTM=DTMF 
IF (A. EQ. DTM) RETURN 
DO 10 1=1,11 
M ( I ) =M ( I ) *A/DTM 
10 CONTINUE 

IF(P5.EQ.0) RETURN 
MW=MW*DTM/A 
RETURN 
END 
C 

SUBROUTINE SOURCE (DX , K , S ,TM , DTM) 
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SOURCE DETERMINES THE VALUE OF THE SOURCE TERM AT A HALF-TIME STEP 

INTEGER P4 
REAL K 

COMMON /AREA6/P3,P4,CF3,TMO,PER,AA,Al,A3,A3 < A4 , A5 

CALCULATION OF HALF TIME 

TMH-TM+ . 5*DTM 
I F ( P4 . EQ . 0 ) GO TO 10 

IF(TMH.GT. AA*PER.AND.TMH. LE . (AA+1 .) *PER)GO TO 10 
AA-AA+1 . 

ADJUSTMENT OF TIME FOR PERIODIC HEATERS 

10 TMP~TMH~AA*PER 

IF (P4 . EQ . 0 ) GO TO 20 
IF (TMP . LE . TMO) GO TO 20 

HEATER OFF 

QF=0 . 

GO TO 30 

HEATER ON 

20 QF=A1*TMP+A2+A3*C0S ( A4*TMP+A5) 

CALCULATION OF SOURCE TERM 

30 S=P3*CF3*DX*QF/K 
RETURN 
END 

SUBROUTINE CON VS (TOLD , T, CMAX, ICV) 

CONVE DETERMINES THE PERCENT DIFFERENCE BETWEEN THE NEW AND OLD 
TEMPERATURES AT A NODE FOR THE GAUSS-SEIDEL METHOD 

DIFF=100.*ABS( (T-TOLD) /T) 

IF (DIFF.GE.CMAX) ICV=1 

RETURN 

END 


C 
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SUBROUTINE WLAYER (TW,TWO,H,HO, T0 ,TO0 ,T1 ,TA2,M, CMAX, ICV, IGS) 

WLAYER CALCULATES THE TEMPERATURES IN THE ICE-WATER LAYER USING 
THE METHOD OF WEAK SOLUTION 

IMPLICIT REAL (K-N) 

COMMON /AREA1/MP , HSMP , HLMP , HMP , CS , CL , TMP , KS , KL 
COMMON /AREA 2/ JW, MW , NW,NA2W 

DIMENSION TW ( 91 ) ,TWO(9i) ,H(91) ,HO(91) ,CNDO(91) 

CALCULATION OF ENTHALPY AND TEMPERATURE AT THE SHIELD-WATER 
INTERFACE 

TOLD=TW ( 1 ) 

I F ( IGS . NE . 1 ) GO TO 9 
CALL PHASE ( 1 , HO ( 1 ) ,HO(2) ,JP0,K02) 

CNDO ( 1 ) =K02* (TWO ( 2 ) -TWO ( 1 ) ) 

K01=K02 

9 CALL PHASE ( 1 , H ( 1 ) , H ( 2) , JP , K2) 

GO TO(10,11,11,12) ,JP 

C ICE 

10 H(1)~(HO(1)+MW*NW*(T0+TO0+ (M-l.) *TWO(l) ) +MW* ( K2*TW ( 2 ) 

1+CNDO ( 1 ) ) )/(l.+MW*(NW*(M+l.)+K2)/CS) 

TW ( 1 ) =H ( 1 ) /CS 
IF (H ( 1 ) . LE . HSMP) GO TO 15 
JP=1 
C M . P . 

11 H ( 1 ) =HO ( 1 ) +MW*NW* (T0- (M+l . ) *MP+TO0+ (M-l . ) *TWO ( 1 ) ) 

1+MW* (K2* (TW ( 2 ) -MP) +CNDO ( 1 ) ) 

TW ( 1 ) =MP 

IF(H(1) .GT. HSMP. AND. H(l) ,LT.HLMP)GO TO 15 
IF (JP.EQ. l.OR. JP.EQ.4)GO TO 15 
C LIQ. WATER 

12 H ( 1 ) = (HO ( 1 ) +MW*NW* (T0— (M+l . ) +TMP+TO0+ (M-l . ) *TWO ( 1 ) ) 

1 +MW* (K2* (TW (2) -TMP) +CNDO (1) ) ) / ( 1 .+MW* (NW* (M+l . ) +K2) /CL) 

TW ( 1 ) =H ( 1 ) /CL+TMP 

IF (H ( 1) . GE . HLMP) GO TO 15 

JP = 4 

GO TO 10 
15 T1 =TW ( 1 ) 

K 1 =K 2 

I F ( ICV . NE . 0 ) GO TO 16 
CALL CONVE(TOLD,TW(l) ,CMAX,ICV) 

JWl-JW-1 

I F ( JW . EQ . 2 ) GO TO 28 


16 
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DO 27 J*2,JW1 
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CALCULATION OF ENTHALPYS AND TEMPERATURES IN THE INTERIOR OF 
THE ICE-WATER LAVER 

TOLD=TW ( J ) 

IW=2 

IF(J.EQ.JWl) IW-3 

IF ( IGS .NE . 1 ) GO TO 19 

CALL PHASE ( IW, HO ( J) ,HO(J+l) ,JP0,K02) 

CNDO (J) =K01* (TWO(J-l) -TWO(J) ) -K02* (TWO ( J) -TWO ( J-f 1 ) ) 
K01=K02 

19 CALL PHASE (IW,H (J) ,H(J+1) j JP,K2) 

GO TO (20 ,21, 21, 22) , JP 

ICE 

20 H(J)s(HO(J)+.5*MW* (K1*TW (J-1)+K2*TW(J+1) -fCNDO ( J ) ) ) 
1/(1.+.5*MW*(K1+K2)/CS) 

TW(J)=H(J)/CS 
IF (H ( J) .LE„HSMP)C,0 TO 25 
JP = 1 
M.P. 

21 H(J)=H0(J)+.5*MW*<K1MTW(J-1)-MP)--K2MMP-TW(J + 1) ) 

1 -fCNDO ( J) ) 

TW ( J) =MP 

IF(H(J) .GT.HSMP.AND.H(J) ,LT.HLMP)GO TO 25 
IF ( JP . EQ . 1 • OR* JP . EQ . 4 ) GO TO 25 
LIQ. WATER 

22 H(J)=(HO(J)f. 5*MW* (Kl* (TW (J-l) -TMP) -K2* (TMP-TW ( J+ 1 ) ) 
1+CNDO ( J ) ) )/(l.+.5*MW*(Kl+K2)/CL) 

TW ( J) =H ( J) /CL+TMP 
IF(H(J) .GE.HLMP)GO TO 25 
JP=4 

GO TO 20 
25 K1=K2 

IF(ICV.NE.0)GO TO 27 

CALL CONVE (TOLD, TW ( J) ,CMAX,ICV) 

27 CONTINUE 

CALCULATION OF ENTHALPY AND TEMPERATURE AT THE OUTER-AMBIENT 
INTERFACE 

28 TOLD-TW(JW) 

IF ( IGS. HE. 1) GO TO 29 

CNDO ( JW) =K01* (TWO ( JW-1 ) -TWO ( JW) ) 

29 K 2 = - 1 . 
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CALL PHASE(2,H(JW) ,H(JW) ,JP,K2) 

GO TO<30, 31,31, 32) ,JP 

* C ICE 

30 H ( JW) « (HO ( JW) +Jfr/*NA2W* ( 2 . *TA2-TWO ( JW) ) +MW* ( K1 *TW ( JW~1 ) 
1+CNDO ( JW) ) ) / ( 1 . +MW* (NA2W+K1) /CS) 

TW(JW)=H(JW)/CS 
IF (H ( JW) . LE . HSMP) GO TO 35 
JP*I 
C M.P. 

31 H ( JW) »HO ( JW) +MW*NA2W # ( 2 . *TA2-MP-TWO (JW) ) +MW*(K1 
1* (TW (JW-1 ) -MP) +CNDO ( JW) ) 

TW(JW)«MP 

IF(H(JW) .GT . HSMP-, AND. H (JW) ,LT.HLMP)GO TO 35 
IF(JP.EQ.1.0R.JP.EQ.4) GO TO 35 
C LIQ. WATER 

3 % H ( JW) - (HO ( JW) +MW*NA2W* ( 2 . *TA2-TMP-TWO ( JW) ) +MW* <K1 
1 * (TW ( JW-1 ) -TMP) +CNDO (JW) ) ) / (1 . +MW* (NA2W+K1) /CL) 

TW ( JW) =H (JW) /CL+TMP 
IF(H(JW) .GE.HLMP)GO TO 35 
JP=4 

GO TO 30 
35 K2=l. 

IF ( ICV .ME . 0) RETURN 
CALL CONVE (TOLD/TW ( JW) ,CMAX,ICV) 

RETURN 
END 

SUBROUTINE PHASE ( IH , HI , H2 ,JJP f K2) 

PHASE DETERMINES THE PHASE AROUND A NODE, AND SETS PHASE 
DEPENDENT CONSTANTS 

IMPLICIT REAL (K-N) 

COMMON /AREA1/MP,HSMP,HLMP,HMP,CS,CL,TMP,KS,KL 
DIMENSION H(2) ,JP(2) ,K(2) 

H (1) «H1 
H ( 2) =H2 
DO 20 1=1,2 

IF (H ( I) . LE . HSMP) JP(I)=1 

IF (H ( I ) .GT. HSMP. AND. H(I) .LT.HMP) JP(I)=2 
IF (H (I) ,GE. HMP. AND. H ( I ) .LT.HLMP) JP(I)=3 
IF (H (I) .GE.HLMP) JP(I)=4 
20 CONTINUE 
JJP=JP(1) 

IF ( K2 . LT . 0 . ) RETURN 
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C 

C SETTING OP PHASE DEPENDENT CONSTANTS 

C 

GO TO ( 25 , 30 # 26) # IH 

25 IF < JP ( 1 ) .NE.2.AND.JP{1) .NE.3)GO TO 30 
H(l) a HMP+.5*(H(l) -HbMP) 

JP (,l) =3 
GO TO 30 

26 IF(JP{2) . NE . 2 . AND . JP ( 2) ,NE.3)GO TO 30 
H ( 2 ) -HMP- . 5* (HLMP-H ( 2) ) 

JP ( 2 ) =2 

30 JPl-JP(l) 

GO TO ( 31 , 31 # 32 # 33) ,JPl 

31 K ( 1 ) ®KS 
GO TO 35 

32 X* (HLMP-H ( 1 ) ) / (HLMP-HSMP) 

K ( 1 ) =2 . *X*KS+ ( 1 . -2 . *X) *KL 
GO TO 35 

33 K ( 1 ) =KL 
35 JP2*JP(2) 

GO TO(41, 42,43.43) #JP2 

41 K ( 2 ) =KS 
GO TO 45 

42 Y= (H ( 2) -HSMP) / (HLMP-HSMP) 

K(2)=2.*Y*KL+(1.-2.*Y) *KS 
GO TC 45 

4 3 K ( 2 ) * 

45 K2- »5*(K (1) +K<2) ) 

RETURN 

END 

/* 

//GO.FT06P001 DD 3YSOUT=A#OUTLIM=9900 
//GO.SYSIN DD * 

II- 5, Pl=l# P2=l # P3=l # P4-1 # P5=l # P6=l # P7 =1 
P8 = 0# P9 = l,Pl0=l,Pll=:l r P12=0,P13=0#P14-0 


JJ- 

15 

10 2 

5 

7 31 

L- 

.087 

.05 .004 

.01 

.012 0.25 

K a 

66.5 

.22 7.6 

.22 

8.7 1.416 

DIF- 

PI- 

HI- 

1.65 

.0087 .138 

.0087 

.15 .0492 

I,Q» 

TMO- 

3# Cl = 
10. ,TMF = 

# C2= 

20. 

25., C3= 

> 

a 

it 

> 

tn 

TIN = 

-4 . , TAl = 

-4 . , TA2= 

-4.# Hl= 

1.# H2= 1000000 

ISI* 

# I SM = 

, DTMI 

. 1 # DTMM 

# DTMF 
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ISP- 250,IFRQ 
IWI= 50,IWM= 
IX= 5 , JX= 

CPS= .5020, KS= 
CPL- .997, KL= 
JW= 31 , JSH= 
***** STANDA R 

/* 


5 , CMAX .005 
150, WI= 1.70, 
7 , TMAX 32.0 
416, DENS 57.4, 
. 32 , DENL 62.4 
I LW= .25 
DE-ICE* DESIGN 



WM= 1.53, WF= 1.30 
MP>= 32., DH= 143.4 


* * + * * 








